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Abstract: We present the complete automation of the universal subtraction for- 
malism proposed by Frixione, Kunszt, and Signer for the computation of any cross 
section at the next-to-leading order in QCD. Given a process, the only ingredient to 
be provided externally is the infrared- and ultraviolet-finite contribution of virtual 
origin. Our implementation, currently restricted to the case of e^e~ collisions, is 
built upon and works in the same way as MadGraph. It is particularly suited to par- 
allel computation, and it can deal with any physical process resulting from a theory 
implemented in MadGraph, thus including the Standard Model as well as Beyond the 
Standard Model theories. We give results for some sample processes that document 
the performances of the implementation, and show in particular how the number of 
subtraction terms has an extremely mild growth with final-state multiplicity. 
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1. Introduction 



The vast majority of new-physics searches at the LHC, and some of those at the Teva- 
tron, are based on final states signatures that feature the presence of a large number 
of jets (typically, from four to ten), generically denoted as multi-jet configurations. 
The fairly difficult problem of giving theoretical predictions for multi-jet processes 
has been completely solved in the past few years, but the solutions are only accurate 
to the leading order (LO) in QCD. All approaches are thus based on the computa- 
tions of tree-level, multi-parton matrix elements, performed by dedicated computer 
programmes (see e.g. ref. [1] for a review and a list of references); the LO accuracy 
implies that the computations can be carried out in four dimensions, and that all 
divergences resulting from integration over the phase space are avoided by means of 
kinematic cuts. Through the computation of the matrix elements, one obtains final 
states composed of quarks, gluons, and other accompanying particles (e.g. VF's or 
Z's). At this point, one may choose to identify each light quark and gluon with a 
jet, thus making use of the local hadron-parton duality, and compare the theoretical 
results to data; we shall call these results matrix- element predictions. A more realis- 
tic approach, which is widely used by experimental collaborations, is that of feeding 
parton-level final states obtained from matrix element computations to event gener- 
ators, which by showering them eventually result in hadron-level final states, where 
any jet may contain several tens of particles; we shall call these showered predictions. 
As is known, showered predictions are only meaningful in the context of a proper 
matching formalism, which avoids the double counting of configurations that can be 
obtained both from the matrix elements and from the showers. Several solutions 
are available for the problem of tree-level matching [2, 3, 4, 5], which involve the 
simultaneous treatment of final states with different multiplicities. 

While fairly successful phenomenologically, the approaches discussed above suf- 
fer from the usual limitations of LO computations. Scale uncertainties, which are 
typically large for processes with high-multiplicity final states, render it unreliable 
the predictions of the absolute values of the cross sections. A common strategy is 
therefore that of fixing the normalization equal to the data for a given jet multi- 
plicity, and of predicting the cross sections for larger and smaller multiplicities. It 
is therefore desirable to extend the accuracy of multi-jet computations to the next- 
to-leading order (NLO) in QCD. It is in fact worth recalling that NLO results, so 
far available only for small-multiplicity jet final states, have been a very important 
ingredient in establishing QCD and the SM as the correct theories of strong and 
electroweak interactions at LEP, SLD, and the Tcvatron, as well as to exclude the 
presence of Beyond the SM signals in the data so far. We also recall that NLO cal- 
culations, apart from resulting in smaller theoretical uncertainties that allow one to 
trust absolute predictions, induce non-trivial effects in particle spectra that cannot 
be obtained in general by a simple rescaling of LO results. 
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The problem of the computation of multi-jet cross sections at the NLO has 
attracted considerable attention in the recent past. There are two main obstacles that 
must be cleared. Firstly, one needs to compute the one-loop corrections. There are 
now several different solutions to this problem, that are reasonably automated, and 
are based on unitarity methods: BlackHat, CutTools, and Rocket (see refs. [6, 7, 8] 
for recent results, and for a more complete list of references). A more traditional 
approach, based on Feynman diagrammatics, is that of GOLEM [9] . It is reasonable 
to expect that most of these codes will attain larger flexibility and speed in the near 
future. Secondly, one must compute the real-emission contributions, and combine 
them with the one-loop ones in order to obtain the physical cross sections. 

The aim of this paper is that of addressing the second of the problems mentioned 
above. We assume that ultraviolet-renormalized one-loop corrections are given, and 
we fully automate the following steps: 

a) computations of real-emission matrix elements and of their local counterterms; 

h) their combinations with the one-loop matrix elements and, if relevant, with 
initial-state coUinear counterterms; and the subsequent definition of finite short- 
distance cross sections; 

c) integration over the phase space of the cross sections obtained in h)-, 

d) output of the results on an event-by-event basis, in the form of a set of four- 
momenta and a weight, that one can use to construct (infrared-safe) observ- 
ables. 

We stress that item h) involves the cancellation of infrared poles as prescribed by the 
KLN theorem. Even if the explicit form of the one-loop contributions is not available, 
this task can be carried out analytically, thanks to the fact that the general structure 
of the poles is known for any process; the residues are proportional to tree-level matrix 
elements that can be computed automatically. 

There are two main ingredients in our work: a programme that computes tree- 
level matrix elements, and a universal formalism for the analytical cancellation of in- 
frared singularities, that allows one to implement in a computer codes short-distance 
cross sections free of any singularities. For the former, we use MadGraph [10, 11]. 
For the latter, we adopt the formalism of Frixione, Kunszt, and Signer, presented 
originally in refs. [12, 13], which we shall refer to as FKS in the rest of this paper. 
Through the implementation of items a)-d) above, we basically achieve an NLO 
version of MadGraph/MadEvent, up to a missing (infrared- and ultraviolet-finite) 
piece of one-loop origin, in which the only operation required from a user is that of 
typing in the process to be computed. We call this programme MadFKS. The current 
implementation of MadFKS, and the results given in this paper, are relevant only to 
the case of e+e~ collisions, which is sufficient to highhght the capabihties of the code. 
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On the other hand, the formahsm is written in full generality, and can therefore be 
immediately applied to any other kind of collisions. We defer to a forthcoming paper 
the complete implementation of the cases where initial-state hadrons are present. 
We stress, however, that what is done here will apply without any modifications (ex- 
cept for the trivial multiplication by the relevant parton density functions) to those 
cases as well. Thus, the only new piece of implementation required by initial-state 
hadrons will concern the emission of a parton from an initial-state parton. We finally 
point out that the results of MadFKS are matrix-element predictions, but also give 
the necessary building blocks to arrive at NLO showered predictions in the context 
of matching techniques such as MC@NLO [14] or POWHEG [15, 16]. 

This paper is organized as follows. In sect. 2 we give a general overview of the 
issues in NLO computations of multi-jet cross sections, and of the defining features 
of the FKS subtraction method. In sect. 3 our notations are introduced. In sect. 4 
we summarize the most relevant formulae of the FKS subtraction for any type of 
collisions, including hadron-hadron ones. In sect. 5 we explain how these formulae 
can be implemented in a computer code. A few refinements which we use in MadFKS 
to optimize the performances of the code are discussed in sect. 6. In sect. 7 we 
present selected results of MadFKS, for the case of e^e~ collisions. Section 8 reports 
our conclusions. Further technical information such as longer formulae, and possible 
variants of the implementation, are given in the appendices. 

2. Computations of multi-jet cross sections at the NLO 

Given a production process, the computation of its NLO corrections in QCD im- 
plies the evaluation of the corresponding one-loop and real-emission corrections (for 
a pedagogical introduction to these concepts, see e.g. sect. 4 of ref. [1]), which are 
eventually added up to get the physical cross section. Although from the principle 
point of view there is no difference between the computation of a small-multiplicity 
and that of a large-multiplicity jet cross section, in practice the situation is entirely 
different. One-loop computations based on straightforward Feynman diagram evalu- 
ation rapidly become too involved to be carried out^, and new techniques have been 
developed recently to bypass this problem (see e.g. refs. [18, 19, 20, 21]). On the 
other hand, the computation of real-emission and Born (i.e., tree-level) diagrams 
poses no problem, and is limited only by CPU. 

There are two issues when summing one-loop and real-emission contributions. 
The first is that of proving (analytically) the cancellation of the infrared singularities 
that arise in the intermediate steps of the computation, and to write the leftovers as 
finite contributions that can be implemented in computer codes. The second is the 

^It is therefore quite remarkable that the authors of ref. [17] have computed the NLO corrections 
to the hadroproduction of ttbb by means of standard techniques. 



4 



actual implementation of these finite terms, which crucially includes an integration 
over the phase space. The former issue has been fully solved [12, 22, 23, 24] in the 
90's, for any jet multiplicity (in other words, the jet multiplicity is simply a param- 
eter in the resulting formulae), in the context of the so-called universal cancellation 
formalisms, which use either the (approximate) slicing method, or the (exact) sub- 
traction method; it is nowadays acknowledged that the slicing method is unsuited 
for describing complicated final states, such as those in multi-jet production. 

The latter of the two issues mentioned above has actually not even been con- 
sidered, for multi-jet cross sections, until recently. There are now several propos- 
als [25, 26, 27] that aim at constructing local counterterms for any given real-emission 
matrix elements (item a) above). What is outlined in items a)-d) is achieved in 
ref. [28] (restricted to massless SM particles), and very recently in ref. [29] (which 
includes the treatment of massive SM particles) . All of these approaches have tackled 
the problem using the universal subtraction method originally introduced in ref. [22] 
(referred to as dipole subtraction henceforth). 

As was discussed in the introduction, in this paper we use the FKS subtraction 
formalism. Although dipole subtraction is the most widely used method for NLO 
matrix-element computations, FKS is the only technique used so far in practice for 
performing NLO showered computations, in the MC@NLO and POWHEG frame- 
works. However, as we shall show in this paper, FKS subtraction has also excellent 
performances in the context of matrix-element computations. While dipole and FKS 
subtraction methods are formally equivalent as shown in ref. [30] , differences between 
the two arise because of the way in which the building blocks for the local countert- 
erms of the real-emission matrix elements are combined (these building blocks corre- 
sponds to the soft, coUinear, and soft-collinear singularities). In the dipole method, 
one emphasizes the role of soft emissions. A structure arises where one sums over 
an "emitted" parton, and two colour partners that exchange the emitted parton; the 
sum over dipoles is thus a sum over three indices. In the FKS method, the emphasis 
is on coUinear emissions, and the structure that emerges is therefore one where there 
is a sum over parton pairs, i.e. a sum over two indices. The consequence of this is 
that, in general, the number of independent subtraction terms is smaller in the FKS 
formalism than in the dipole formalism. 

One of the key features of FKS subtraction is that, in the integration of the real- 
emission matrix elements, one effectively defines partonic processes with at most 
one soft and one coUinear singularities, which results in a much simpler subtraction 
structure than that of the original matrix elements. These processes are furthermore 
fully independent from each other, and can therefore be computed separately. This is 
a genuine parallelization, more efficient than that of computing the same production 
process several times with different seeds for random number generation, especially 
in view of the use of adaptive-integration routines (the latter type of parallelization 
can obviously still be set up in the context of FKS subtraction, if need be). 
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The original FKS papers [12, 13] addressed the case of processes with final states 
composed only of massless quarks and gluons. The extension to generic processes, 
with the presence of both strongly-interacting massive particles and of non strongly- 
interacting particles, is almost trivial. Several of these cases have in fact been im- 
plemented in MC@NLO using FKS subtraction (see e.g. ref. [31]). Here, we shall 
present a complete summary of the relevant formulae, which can also be applied to 
computing cross sections to NLO accuracy for the production of non-SM particles. 

In summary, we advocate the use of FKS subtraction for matrix-element predic- 
tions for a variety of reasons: structure identical to that of collinear emissions, thus 
lending itself naturally to matching with parton shower Monte Carlos not based on 
colour dipoles; ease of importance sampling; small number of subtraction terms, and 
the modest growth of their number with the multiplicity, with beneficial effects on 
numerical stability; organization of the calculation in a way which is parallel in na- 
ture; and the availability of all results necessary for the treatment of fully-polarized 
processes, which allows one to performing sums over helicity states with Monte Carlo 
methods. 

In the rest of this paper, all formulae for NLO cross sections will make use of 
FKS subtraction. 

3. Notation 

3.1 Partonic processes 

For any given partonic process, we shall denote by n the number of final-state par- 
ticles at the Born level. The contributions to the NLO partonic cross sections will 
therefore either have a2— >nora2— >(n+l) kinematics, which we shall call Born 
(or n-body) or real-emission (or (n -|- l)-body) kinematics respectively^. They will 
be denoted as follows: 

^1 + ^2 — * /c3 + ■ ■ ■ + kn+2 Born kinematics ; (3.1) 
ki -\- /c2 — ^ ^3 + • — I- ^n+3 real emission kinematics . (3.2) 

We shall write the corresponding phase spaces as 

d(t)n , d4>„+i (3.3) 

respectively. We shall denote by n0 the number of final-state particles which are not 
strongly interacting (e.g. leptons); uh will denote the number of massive, strongly- 
interacting particles (e.g. heavy quarks). Depending on the process, some or all of 
the niD + riH particles will not belong to the Standard Model, n- and (n -|- l)-body 
processes will be characterized by the same values of n0 and of uh'i the identities of 

^The decay of a particle can also be described by using eqs. (3.1) and (3.2). 
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the corresponding particles will also be the same in the two classes of processes. On 
the other hand, the number of light quarks and gluons in the (n + l)-body processes 
(denoted by n^^) will be equal to that in n-body processes (denoted by n^L ^), plus 
one. We have therefore: 

nf = ) + 1 , (3.4) 
n = Ti-p +nH + n$, (3.5) 
n + 1 — n'P + + n0 . (3.6) 

Ifl<A;<n + 2orl<A;<n + 3isan index which runs over all particles, we shall 
adopt the following labeling scheme: 

1 < A; < 2 =^ 

initial state ; (3.7) 
2><k< nf^^ +2 =^ 

masslcss quarks and gluons; (3-8) 

strongly-interacting massive particles ; (3-9) 
u^l'"^^ + + 3 < A; < ^^^^ + + n0 + 2 =^ 

non strongly-interacting particles. (3.10) 

In order to deal with the various types of colliding particles, it will also be useful to 
introduce the quantity n/, which can take the following values 

ni — 1 <^=^ HH collisions, (3-11) 
nj — 2 <^=^ eH collisions , (3-12) 
ni — 3 <^=^ ee collisions , (3.13) 

where e and H generically denote a non-hadronic and a hadronic particle respectively; 
rii is therefore the smallest value that a label can assume when running over light 
quarks and gluons (without loss of generality, we have conventionally given label 1 
to the non-hadronic particle in an eH collision). 

The identity of particle k will be denoted by 2^; it is also convenient to denote 
the anti-particle of by X^, although X^ can be used as well. A process will be 
unambiguously identified by giving the list of the identities of the particles involved, 
and we shall denote it by 

r = (Xi, . . . X„+2) or r = (Xi, . . . X^+g) . (3.14) 

The sets of all n- and {n + l)-body processes will be denoted by TZn and by 7^„+i 
respectively^. In the following, we shall make use of the fact that a given n-body 



^If two processes in the form of eq. (3.14) are related by a permutation over the identities of 
final-state particles, only one of them will be included in TZn or TZn+i ■ 
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process can be related to at least one (n+l)-body process. We therefore introduce the 

( R) 

following notations, where i is a final-state light quark or a gluon (3 < i < u'l' + 2), 
and j is a strongly-interacting particle {nj < j <ni^ + uh + 2): 

r = (Xi, . . .Xj, . . .Tj, . . .Tn+s) £ T^n+i (3.15) 

=^ r^,r^®^'\enn if r\r^®^'^ exist. (3.16) 

We have defined 

= {Il,...Xi:...Ij:...In+3) , (3.17) 

H®^'^ = (Ti, . . .Xi, . ..Im, ■ ■ -^n+a) . (3.18) 

Thus, process r-^ in eq. (3.17) is constructed by simply removing parton i from the 
original list. It is clear that such an operation results in a physical process (i.e., 
which exists) only when i is a gluon. Process r^®*'^ in eq. (3.18) is constructed by 
removing parton i from the original list, and by replacing particle j with one whose 
identity is that of the parton entering the (Jj^j, Jj,2j) QCD vertex (thus, ii li — g 
and Ij — g, then Ij^^ — g; ii Ij, — g and Ij — q, then Ij^i — q, and so forth). If this 
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QCD vertex does not exist (as e.g. in the case Ti = u and Xj = d), process r- 
non physical, i.e. it does not exist^. The definitions of r-^ and of r-'®*'^ are obviously 
motivated by the relevant soft and coUinear limits respectively, but will be used in 
several different contexts^. 

3.2 Matrix elements 

NLO cross sections receive contributions from both tree-level and one-loop ampli- 
tudes, which we shall denote as follows: 

^("+1,0) _^ real- emission tree amplitude (r G TZn+i) , (3.19) 
^(",0)^^^ — ^ Born tree amplitude (r G 7^„) , (3.20) 
^{'».i)(^) — , one-loop amplitude (r e 7^„) . (3.21) 

The amplitudes above are typically relevant to a given spin and colour configuration; 
they include all the coupling constant factors. Starting from amplitudes, one needs 
to construct several quantities which will enter short-distance cross sections. We 
denote them as follows: 

^'""'°'w-s:;(i5^EM'"^"'Mr. (3.22, 

colour 
spin 



^If j = 1 or j = 2, one needs to use 1^ instead of Ti when constructing Tj^i; in the following, we 
may therefore use the notation Ti^i or 22ei- 

^Note that eq. (3.17) is a particular case of eq. (3.18), since the two coincide when Xj = g. 
However, it will turn out to be convenient to have the two different notations. 
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colour 
spin 

•^^"'"'^"^ = ~ L4mX2) ^ '4(-°)(r)g(X.)-g(X0^("'°Hr)*, (3.24) 

colour 

spin 

= S::(2oW ^ {-^'-''Vj^'-'Hr)*} . (3.25) 

colour 
spin 

In these equations, s = {ki + k2)'^ = 2ki-k2, and uj{I) is the product of spin and colour 
degrees of freedom for particle X. In 4 — 2e dimensions, we have uj{q) = 2Nf. = 6 and 
oj{g) = 2(1 — e)Dj^ = 16(1 — e). These average factors serve the sole purpose of fully 
specifying the divergent part of the one-loop contribution in the CDR scheme, which 
is not used in numerical calculations (see app. B for more details). We stress that the 
formulae which give the physical cross sections in FKS, and that are implemented in 
computer codes, are finite in four dimensions. Note that symmetry factors accounting 
for the presence of identical particles in the final state are not included in the matrix 
elements defined above; they will be inserted later in the expressions for the short- 
distance cross sections. As the notation suggests, r e TZn+i in eq. (3.22), while 
r e TZn in eqs. (3.23)-(3.25). The quantities defined in eqs. (3.22)-(3.25) are sufficient 
to perform all of our computations, with the exception of an azimuthal contribution 
to the coUinear limit of which requires a reduced matrix element that we 

shall denote by The precise form of this quantity is irrelevant in this section; 

it will be used in sect. 5.3, and its definition will be given in app. D. 

The colour operators Q{T) that enter the definition of the so-called colour-finked 
Bern's, eq. (3.24), give a representation of the colour algebra, and can be defined as 
follows: 

Qi^) = , {-^"^iLi . {^"}Li , ^ e 3, 3, 8, (3.26) 

with and the SU(3) generators in the fundamental and adjoint representations 
respectively. Clearly, eq. (3.26) can be extended to higher-dimensional representa- 
tions if need be. It is trivial to show that 

g(Xi) • Q(l2) = Q(X2) • Q(Ii) , (3.27) 

g(x) • g(x) = Q\i) = c(x)7, (3.28) 

where 

C{I)^Cf^^^ for Xe3,3, (3.29) 

C(X) = = A^c for X e 8 . (3.30) 



9 



Using the colour-singlet condition 



J2 Qi^k) = (3.31) 

k=ni 

and eqs. (3.27)-(3.30), it is easy to prove that 

Mtl''^ = Mjr^ , (3.32) 



J2 = 2C(T0A^^"'°^ , (3.33) 

M^f = -C(Jfc)A<('^'0) , (3.34) 



k^l 
k=ni 



with 

ni<k,l< njf ^ +nH + 2. (3.35) 

Note that is an interference term related to the exchange of a soft gluon 

between particles k and I, which appears in the soft limit of the real matrix elements. 
This term has a coefficient proportional to an eikonal factor, which therefore vanishes 
identically when k — I and k is massless (see app. A). 

The final formulae of the FKS subtraction method, to be given below, will only 
deal with non-divergent quantities. The interested reader can find in the original 
paper (ref. [12]) the proof of the cancellation of the infrared and collinear singularities 
that arise in the intermediate steps of the computation. As in all computations at 
the NLO, the finitcncss of the partonic short-distance cross sections is in general 
a consequence of imposing kinematic cuts on final-state particles. Without loss of 
generahty, we can always assume these cuts to be equivalent to the request of having 
either or -|- 1 jets in the final state, the jets being reconstructed with an 
arbitrary algorithm. This condition will be symbolically denoted by 

(B) 

. (3.36) 

(S) 

The condition in J"^ is sufficient to prevent the appearance of phase-space singu- 
larities^ in n-body quantities such as Al^"'^'', and TW^"'^). On the other hand, 
the {n -\- l)-body tree-level matrix elements will still diverge in some regions of the 
phase space; these divergences are subtracted (in all subtraction methods) by means 
of suitable counterterms. In order to classify these divergences and to eventually 



°It is understood that the condition in J"^' also prevents aU divergences due to final-state 
leptons and photons. This is necessary if the subtraction of the corresponding QED singularities is 
not carried out. 
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subtract them, for any given process r G TZn+i we introduce the following set of 
ordered pairs (called the set of FKS pairs) 

VFKs{r) = [it,j) I 3<i<n^P +nH + 2,ni<j <nf^ +nH + 2,iy^j , 

_yy^(-+i,o)(^ynf ) ^ ^ ^0 ^ Q ^0 ^ or || a!,}. (3.37) 

In words, a pair of particles belongs to the set of FKS pairs if they induce soft or 
coUinear singularities (or both) in the {n + l)-body matrix elements, even in the 
presence of jet cuts. The first element of the pair will be called the FKS parton, 
and the second element will be its sister^. Conversely, by definition Pfks takes into 
account all phase-space singularities of the (n + l)-body tree-level matrix elements 
after jet cuts. We note that the condition kj ^ is simply not relevant if j = 1, 2, 
since in such a case the momentum kj is fixed. In the computation of an NLO cross 
section according to the FKS formalism, each pair belonging to T^fks will correspond 
to a set of subtractions of soft and coUinear singularities which, when combined with 
the real-emission matrix element, will result into a finite contribution to physical 
observables. Each of these contributions is separately finite, and can therefore be 
computed independently from the others, and integrated with numerical methods. 

It is immediately obvious that some of the pairs with i,j < n"^^ -\- uh + 2 
will not belong to "Pfks foi' any physical process. This is the case e.g. when i and 
j are both massive and strongly interacting; or when i and j arc a quark and an 
antiquark of different flavours; or when one of the members of the pair is massive 
and strongly interacting, and the other member is not a gluon. In general, however, 
for large multiplicities the number of elements in Vyks will scale approximately as 
n"^\n"^^ + t^h)- On the other hand, depending on the identities of the particles in 
r, some of the contributions due to the various FKS pairs will actually be identical, 
because of the symmetry properties of matrix elements and phase spaces. This 
implies a drastic reduction of the number of independent terms actually needed in 
the computation. We shall return to this issue in sect. 6.1. 



^Although in the definition of 'Pfks the role of the FKS partons and of their sisters is symmetric 
if they both belong to the final state, this will not be the case when the subtraction of singularities 
will be performed: hence, the convenience of distinguishing them. 
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4. Cross sections 



We write an NLO cross section in the collision of two particles Pi and P2 with 
momenta Ki and K2 using the factorization theorem 

dap^P^iKi, K2) = 

X (da^^^+^Xm; XiKi, x^K^) + da^^+^\rR, x^K^, 
+ J2 J dxidx2fi:'^'\xi)fi;^'\x2)da^"\rB;xiKi,X2K2). (4.1) 

As the notation suggests, we have (see eqs. (3.1) and (3.2)) 

h = XiKi , k2 = X2K2 . (4.2) 

If Pa is a hadron, /j^"'' will be a parton density function (PDF). Otherwise, it may 
simply be equal to 5{1—Xa) (e.g. when describing an electron beam at a fixed energy), 
or to a more complicated function that effectively describes the energy loss when the 
incoming particle Pa turns into the particle that enters the hard reaction (e.g., 
f^j' may be the Weizsacker-WiUiams function with P^ — e and 2q, = 7). These 
details are in any case irrelevant in what follows, where we shall deal only with the 
short-distance partonic cross sections da^"), da^''+^\ and da^''+^'^] if necessary, we 
shall refer to functions /j^" gencrically as PDFs. 

In this paper, we shall denote by the common value of the factorization and 
renormalization scales, ^ = fip = fip, which simplifies the writing of the formulae. 
At the end, one is always able to recover the separate dependence upon /ip and /ir, 
by exploiting the renormalization group invariance with respect to these two scales. 
All the relevant formulae are given in app. C. 

The short-distance cross sections da^'^\ d(7^"'+^\ and da'^'^~^^\ have an n-body, an 
{n + l)-body, and a degenerate (n + l)-body kinematics respectively. They will be 
discussed in the next three subsections in turn. 

4.1 n-body contributions 

We decompose the n-body cross section into four terms: 

da^^^ = da^^'") + d(7(^'") + da^^'") + da^^-") . (4.3) 
The first term on the r.h.s. of eq. (4.3) is the Born contribution (r € TZn) 

(B) 

Al [r) 
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with M{r) the symmetry factor that takes into account the presence of identical 
particles in the final state. The remaining three terms on the r.h.s. of eq. (4.3) 
have, roughly speaking, a coUinear, soft, and one- loop origin respectively; they are 
separately finite, and we understand the kinematics given in eq. (3.1). We have: 



where 



(B) 

|Q(.)A.(-)(.)^, 



(4.5) 



Q(r) = - log ^ 7(Xi) + 2C(Xi) \ogUt + l{:^2) + ^CiX^) \ogUt 



fc=3 



s5, 



y (Xfe) - log S ( liXk) - 2C{Ik) log 



2Q^ 



+ 2C{Xk) log 



2^ 



cut 



27(2-fe)log 



2-5'fc 

7i 



(4.6) 



Note that, in this equation, only light quarks and gluons are involved. We have 
s — 2ki-k2, and Ek is the energy of parton k in the cm. frame of the incoming 
particles Ii and I2 (i.e., the partonic cm. frame in the case of hadronic collisions). 
The Casimir's C(T) have been defined in eqs. (3.29) and (3.30); the other colour 
factors are: 



1(9) 
7(g) 



6 
3 

2 







67 


27r2 


9 


3 


13 


27r2 




3 



Ca - '^TpNf 
9 



(4.7) 
(4.8) 
(4.9) 

(4.10) 



Equation (4.6) contains two free parameters, ^cut and 5o- The same parameters will 
be used in the subtraction of the soft and final-state coUinear divergences, respec- 
tively, that affect the real-emission contribution, as it will be described below. Thus, 
although eq. (4.5) does depend on ^cut and 5o, the physical cross section does not. 
This constitutes a powerful check of the correctness of the implementation of the 
formalism; we shall return to this point in sect. 7. In eq. (4.6) we have denoted by 
Q the Ellis-Sexton scale [32], which is an arbitrary mass scale that one may use to 
write the one-loop results; more details on its role will be given in app. B and in 
app. C. Here we just stress that the NLO cross section is exactly independent of 
Q. Finally, we remind the reader that /i — /ip — /ir. Note that the contribution in 
round brackets in the first line on the r.h.s. of eq. (4.6) is non trivial only if at least 



13 



one of the incoming particles is a hadron. Still, the form in eq. (4.6) applies to any 
kind of collisions, provided that one defines 

C{I) = 7(T) = if T is a colour singlet. (4.11) 

The n-body contribution of soft origin reads: 

n'-f'>+nH+2 n'lhnH+2 (B) 

da^'-\r) = '^ E E 4r'"^'^>tirW^^0n. (4.12) 

k=nj l=k ^ ' 

The quantities f;^™*^'™') are the finite parts of the integrated cikonal factors; their 
precise definitions and explicit forms are given in app. A. We have denoted by 

mi = ^|k\ (4.13) 

the mass of particle /. We point out that the contributions with k = I and k < n)^ +2 
to eq. (4.12) are identically equal to zero, owing to eq. (A. 4). 
Finally, the n-body contribution of one-loop origin is: 

(B) 

d<.(^.")(.)^|^V&i;(r)^#„. (4.14) 

The quantity V^pj^ is the non-divergent part of M.^'^'^^ introduced in eq. (3.25). Its 
definition is therefore unambiguous only in a definite scheme, and the one to be used 
in eq. (4.14) is the Conventional Dimensional Regularization (CDR) scheme. The 
extraction of Vp^j]J from a complete one-loop computation, and its relation to the 
popular Dimensional Reduction scheme, is discussed in app. B. 

4.2 {n + l)-body contributions 

The short distance cross sections dcr*^"^^^ in eq. (4.1) are due the contributions of the 
real-emission matrix elements ^("+^'°), with their phase-space singularities suitably 
subtracted. As is known, the integration of the subtracted matrix elements is the 
most involved from the numerical viewpoint, and the difficulty increases with the 
number of external legs, because of the proliferation of singularities. In the FKS 
formalism, this problem is simplified by effectively partitioning the phase space, in 
such a way that in each of the kinematic regions resulting from the partition at 
most one soft and one coUinear singularity are present. This partition is achieved by 
introducing a set of positive-definite functions^ 

Sijir), (^,J) ePFKs(r), (4.15) 



*In the original papers, refs. [12. 13], the S functions were constructed using Heaviside 8's. 
Smoother S functions, resulting in an improved numerical behaviour, were introduced in ref. [31]. 
With © functions, one achieves an exact phase-space partition, in which the various regions do not 
overlap, while with the present form the partition regions can overlap. 
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where the argument r G TZn+i stresses the fact that different choices for S can be 
made for different processes; this property is useful to reduce the computing time. 
In the following, we shall often understand the argument r in the S functions. 

There is ample freedom in the definition of the S functions, but the following 
constraint must be obeyed: 

= (4.16) 

(j,j)e'PFKS 

Furthermore, these functions must have the following behaviours in the kinematic 
configurations associated with the singularities of 

lim Sij = hiji ) if rrii ^ rrij ^ , (4.17) 

k,\\kj + 

lim Sij — Cij if Xi = g , with 

fe9— ►0 

< Cij < 1 and ^ = 1 . (4.18) 

i 

(j,i)e'PFKs 

lim Sij = y{k,l}y^ {ij} with (A;, I) e Pfks and ruk ^ mi ^ 0, (4.19) 
lim Sij = VA; with Tk = g and 3/ with {k, I) e Vfks or {l,k) e Vfks ■ 

(4.20) 



In words, Sij goes to zero in all regions of phase space where the real emission matrix 
elements diverge, except if this involves particle i being soft, or particles i and j 
being coUinear. Note that eqs. (4.17) and/or (4.19) need not be satisfied if one of 
the partons in the relevant FKS pair is massive; in such cases, the limits can assume 
arbitrary values. Also, eq. (4.20) implies in particular that the limit of Sij is zero 
when j is soft; thus, although according to eq. (3.37) there may be soft singularities 
associated with FKS sisters, they are damped by using Sij. The functions hij{z) 
introduced in eq. (4.17) are defined in < 2; < 1, and have the following properties: 

hij{z) = 1 if n/ < J < 2, (4.21) 

hij{z) = h{z) if 3<j< nf^ + n^/ + 2 , (4.22) 

with h{z) a positive-definite function such that 

lim/i(^) = l, lim/i(^)-0, h{z) + h{l - z) = I . (4.23) 

The physical cross section is independent of the choice of h{z)] we shall return to 
this point in sect. 7. As stressed above, the particular form of the S functions is not 
important in what follows; we shall give an explicit construction in sect. 5.2. As a 
final remark, we point out that eqs. (4.21) and (4.22) imply that the notation hij is 



15 



redundant here; it will however become useful in the context of the optimization we 
shall carry out in sect. 6. 

Using eq. (4.16), one proceeds by writing 

^{n+i,o)(^^^^ 5y(r)M("+^'°)(r) . (4.24) 

Thanks to eqs. (4.19) and (4.20), each term on the r.h.s. of eq. (4.24) will only be 
singular when particle i is soft, and/or particles i and j are coUinear^. Therefore, each 
term can actually be regarded as describing a production process with the simplest 
possible structure of phase-space divergences. Furthermore, these contributions are 
fully independent from each other. 

We shall now give the expressions of the subtracted real-emission cross sections. 
In order to do so, we work in the cm. frame of the incoming particles: 

h = ^(1, 0, 0, 1) , ^2 = ^(1, 0, 0, -1) . (4.25) 
In this frame, for each pair (i, j) e Pfks we introduce the variables and i/ij, where 

Ei^^^^, (4.26) 
Vij . (4.27) 

In other words, is the rcscaled energy of the FKS parton, and yij is the cosine of 
the angle between the FKS parton and its sister. Using these variables, the soft and 
coUinear singularities of SijAd'^"^^'^^ correspond to = and to i/ij = 1 respectively 
(note, however, that the matrix element may not be singular in either or both of 
these limits, depending on the identities of particles i and j). We have (see ref. [12]) 

da(^+'\r)= J2 da^r'\r), (4.28) 

(i,i)€PFKS 

where 

(4.29) 

The variable ipi parametrizes the azimuthal direction of the FKS parton, but its 
precise definition is irrelevant here. The quantity is related to the n-body phase 
space, and is implicitly defined as follows: 

d(j)n+i = ^id^idyijdipid^'^ . (4.30) 
^We remind the reader that multiple soft and collinear singularities, i.e. configurations that 

(B) 

would contribute to NNLO, are cut by J"^ . 
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The actual form for d(f)l{ depends on the specific parametrization adopted for d<pn+i, 
but the following equations must always hold: 

lim d^'l = 7T^d(l>n if = , (4.31) 
lim d(^^ = -r;^d4>n if = m,- = . (4.32) 
The distributions entering eq. (4.29) are defined as follows^°: 



1 



1 -/-l 1 Vij 



where 



1 \ 

1-- J]mJ . (4.35) 



>,fe=3 

In eq. (4.33) and (4.34) ^cut and S are free parameters, that can be chosen in the 
ranges 

< < Uax , 0<6<2. (4.36) 

It would be possible to choose different values of C,cut and 6 for each of the cross 
sections da^J^^^ contributing to eq. (4.28). In practice, we shall choose the same 
value of ^cut for all {i, j) pairs, the same value of 6 for all pairs (i, j) with j > 3 (which 
we shall denote by Sq), and the same value of 5 for all pairs with j <2 (which 
we shall denote by Si). The parameters ^cut, So, and Sj are therefore associated with 
soft, final-state coUinear, and initial-state coUinear singularities respectively. As was 
discussed in sect. 4.1, the physical cross section on the l.h.s. of eq. (4.1) is strictly 
independent of ^cut, So, and Sj, while this is not the case for the short-distance cross 
sections rfa^"), da'-''+^\ and da^^'+^l 

For the reader unfamiliar with the use of plus distributions in the context of 
cross section calculations, it is useful to explicitly show the structure of eq. (4.29). 
We start by rewriting such an equation with the following shorthand notation: 

d4^'^ = (t^) , Vij)d^idVii ■ (4-37) 

We can now expand the plus distributions using their definitions, eqs. (4.33) and (4.34): 

- / dii 1 dyij— Eij(^i, yij) - l)Q(yij - I + S) 

Jo J-l SiU ~ Vij) L 

- E,, (o, VijMUt - Ci) + ^iM ^MUt - ^iMvij - 1 + <^)] ■ 

(4.38) 



^°If rrii ^ 0, the lower limit of the integration range is larger than zero. This case is however 
trivial, and we shall return to it at the end of this section. 
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There are four terms in the integrand in eq. (4.38). We shaU caU them "event" 
{'Eij[$^i,yij)), "colhnear counterevent" (Sjj(,^j, 1)), "soft counterevent" (Sjj(0, t/jj)), 
and "soft-colhnear counterevent" (Sjj(0, 1)). In general, they may have four different 
kinematic configurations, but one can actually reduce them to two, as discussed in 
sect. 5.1. Note that, in the case of a massive sister {mj ^ 0), the colhnear and soft- 
collinear counterevents vanish identically, because of the damping factor (1 — y^j) in 
eq. (4.29). 

Finally, consider the case ^ 0. According to eq. (3.37), (i,j) e PrKS only 
if Xj = g. The only possible singularity arising from the FKS pair is the 

soft one kj 0; this is however damped by Sij. It follows that all counterevents in 
eq. (4.38) vanish, and the distributions in eq. (4.37) coincide with ordinary functions. 
Therefore, using eq. (4.30), one simply gets 

(S) 

dal;^'\r) = M(-+i'0)(r)5,,(r)-^ , m, ^ and Ij ^ g . (4.39) 

Since there are no subtractions involved, one is not bound to use the variables ^j, i/ij, 
and (fii in the parametrization of d(f)n+i in this equation. On the other hand, there is 
nothing that prevents one from doing so; in such a case, note that the smallest value 
can assume is 2mj/y^, and we understand in eq. (4.38) that T^ij^^i.yij) vanish 
identically if < 2mj/ y^. The fact that when ^ one does not need to subtract 
any singularities suggests that from the viewpoint of the efficiency of the calculation 
it is not convenient to define the quantity d(Tlj~^^^ in eq. (4.39) as an independent 
contribution; we shall discuss this issue in sect. 6. 

4.3 Degenerate {n + l)-body contributions 

These contributions, denoted by da'^'^'^^^ in eq. (4.1), are the finite remainders left by 
the subtraction of initial-state colhnear counterterms. We write 

da(-+^)(r)= Yl Mr\)+ E M2'-'\r). (4.40) 

(j,l)ePFKS (i,2)€PFKS 

The condition G "Pfks (with j = 1,2) in eq. (4.40) implies that the real- 

emission matrix clement is singular when ki \\ kj. This typically means that, if Pj 
is not a hadron, the corresponding contribution da^J^^^ need not be included. A 
notable exception is Ij = g and 2j = 7, when the matrix element has a (QED) 
singularity, which is cancelled by the counterterm due to the inhomogeneous part 
of the (Altarelli-Parisi) evolution equation for the hadronic photon. In any case, by 
writing (i, 1) G T'fks ^-nd {i, 2) G T'fks in eq. (4.40) we can use this equation for any 



18 



kind of collisions. As shown in ref . [12] , one obtains 



dali'^^\r] ki, k2) 



27r i 



(S) 



s6i 



■6) 



■6) 



log 6 



(4.41) 



-i^X23.X2(l-Ci) 



ki,{l-Ci)k2)dCi, (4.42) 



where, in analogy with eq. (4.33), we have introduced the distribution 

'^^6/(6) (^) =^^""^''c(/(6)-/(o)e(u-eo)^- (4.43) 



We have also defined the functions 



Pa,iz, e) = (1 - z)P,,{z, e) = Pi?(z) + eP^^iz) + Die') , 



,75(1), 



(4.44) 



i.e. the unregularized Altarelli-Parisi splitting functions times a {1 — z) factor that 
damps the z ^ 1 soft singularity. The explicit forms of Pj^'^^ and P^l^ are given 
in app. D. The functions Kab are related to the choice of the PDF scheme; for all 
practical purposes they are trivial, since in the MS scheme we have Kab — 0. We 
point out that in the special case alluded to at the beginning of this section, i.e. li — q 
and Ij — 7, the relevant da^^^^^ must be proportional to aem- We can therefore still 
use eqs. (4.41) and (4.42), provided that we include a factor aem/ots in the kernels 
Pq^ and Kq.^ (see ref. [13]). 



Consider da^^l^^^ in eq. (4.41) (the case of da'^l^^^' is fully analogous). The fact 
that (i, 1) G "Pfks implies that r^®*'^ exists and belongs to 7^„, and that Ti^f^ is either 
a light quark or a gluon. The notation embodies the fact that the matrix element 
and phase space must be computed with 2 — > n kinematic configurations, where 
incoming partons have momenta 



k[^{l-ii)k^, k'2^k2. 



(4.45) 



We conclude by stressing that the symmetry factor in eqs. (4.41) and (4.42) 
is Af{r), and not J\f(r^®'^'^) or J\f{r'^®'^'^). This can be easily understood as follows. 
Suppose that there are p identical particles in the final state of r, with labels ii, . . .i 



pi 
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jV(r) will thus contain a factor p\. On the other hand, we shall also have da^^^^^ — 
d^ut^\ for j — 1, 2, and I < a, /3 < p. Hence, 

t<T'=J'<""- (4-46) 

a=l 

Apart from allowing one to reduce the computing time by a factor of p, this equation 
also shows that symmetry factors j\/'(r^®'i'^i) and jV^r^®'!'^^) will emerge naturally, 
since 

(4.47) 



Afir) p\Jf{r) {p-l)\Jf{r) A^(ri®^i'V) Af{r^®ii,\i) 
This equation obviously holds also for p = 1. 



5. Implementation 

The most involved part in the computation of an NLO cross section is the integration 
of the subtracted real-emission (i.e., (n-l-l)-body) matrix elements. In this section we 
shall therefore mainly present technical details on the ingredients that will enter such 
an integration in our implementation. In particular, we shall discuss the kinematics 
relevant to the generation of count ere vents, give the formulae relevant to the evalu- 
ation of the damped matrix elements in the singular limits, and present an explicit 
construction of the S functions. We stress that the latter is fully general, and easy 
to automate, but nothing prevents one from using alternative forms, provided they 
fulfill the conditions given in eqs. (4.16)-(4.20). Although the practical applications 
in this paper are limited to the case of e'*"e~ collisions, the formulae we give in the 
following can be applied without modifications to any kind of incoming particles. 

5.1 Kinematics 

The integration of the subtracted real-emission matrix elements over the (n -|- 1)- 
body phase space entails the generation of 3n — 1 random numbers, of which three 
correspond to the integration variables yij, and (fi, that have been singled out since 
they are directly associated with the FKS parton. We denote the remaining 3n — 4 
random numbers by Xa] we can write a generic 2 — > n -|- 1 kinematic configuration 
as follows: 

|- -] n+3 

As imphed by eq. (4.38), this configuration will be that of the event. Starting from it, 
we construct the configurations corresponding to the soft, coUinear, and soft-coUinear 
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counterevents 



,11. 




n+3 



(5.2) 



(5.3) 



(5.4) 



with (fii and in eqs. (5.2)-(5.4) the same as those used in the generation of the 
kinematics of the event, eq. (5.1). For a given choice of the 3n — 1 independent 
variables, the momenta ki depend on the phase-space paramctrization adopted. We 
have considered three different such parametrizations, in order to test the convergence 
behaviour and numerical stability of the numerical integration of short- distance cross 
sections. All obey the following equations: 



These properties trivially follow from the observation that the kinematic configu- 
rations of the counterevents are degenerate, and effectively correspond to a 2 — >• n 
configuration. This configuration, whose momenta we denote [14, 16] by k, can 
simply be defined with a relabeling of the soft counterevent kinematics: 



with a denoting the relabeling. 

We stress that, although eqs. (5.5) and (5.6) need not be fulfilled by any phase- 
space parametrization, and are not mandatory for the implementation of FKS sub- 
traction, they are quite intuitive from the physical point of view, and help improve 
the numerical stability of the results for differential distributions, while also reducing 
the computing time. 

We finally point out that the results of this section apply to the case j > 3, i.e. to 
final-state emissions. The situation is more involved in the case of j = 1,2, i.e. of 
initial-state emissions, but similar simplifications are possible (see e.g. app. A. 4 of 
ref. [14]). We postpone a discussion on this issue to a forthcoming publication. 

5.2 Choice of S functions 

As already mentioned in sect. 4.2, in ref. [31] a smooth form of the S function was 
adopted, which was seen to improve the numerical integration with respect to the 
form of refs. [12, 13], where Heaviside O's had been used. In ref. [31], the S functions 

^^In the case in which particle j is massive, there is no need to construct coUinear and soft-colUnear 
counterevents. If i is massive, all counterevents are identically equal to zero. 



(5.5) 
(5.6) 




(5.7) 
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were constructed using invariants k^-ki, but given that FKS subtraction uses energy 
and angles as integration variables, in a further generalization (see e.g. ref. [16]) one 
can also define: ^ 

^ij = hij{zij), e Vfks , (5.8) 

where the functions hij have been introduced in eqs. (4.21)-(4.23), and 

Ei El 

~ P _1_ F ' ~ 77. I 77. ~ -'- ~ % ' (5-9) 



V ^hui{zui), (5.10) 

{k,l)ePFKS 

dki^ f^y'(l-/5,Acos^,0'' , (5.11) 



Pk-Jl-^. (5.12) 



The parameters and 65 introduced in eq. (5.11) are real, positive, and arbitrary. 
The fact that the physical cross section is independent of their choices constitutes 
another check of the correctness of the implementation of FKS subtraction. We shall 
return to this point in sect. 7. As shown in eqs. (4.21) and (4.22), the functions hki 
are given in terms of a function h{z), that we have chosen as follows: 

(I — zV^h 

h{z) = „ \ ,/ , an^l. (5.13) 

In practice, the specific form of h{z) has only a modest impact on numerical compu- 
tations; we shall briefiy comment on this point in sect. 7. The role of h^i in eq. (5.10) 
is simply that of avoiding to count twice the contribution of l/d^i when both (A;, I) 
and (Z, k) belong to T'fks- In f^-ct, in this case one can immediately obtain 

-^hki{zki) + -^hik{zik) = (hki{zki) + hki{l - Zki)) = -j"- , (5.14) 

dkl dlk dkl ^ ^ dkl 

having used eqs. (4.23) and (5.9). Alternatively, one could have dropped the hki 
factor in eq. (5.10), and summed over unordered pairs. 

We point out that the limits of the <S functions given in eqs. (4.17)-(4.20) are 

related to the vanishing of dij in the coUinear (6ij with rrii = rrij = 0) and soft 
{Ei 0) limits. When parton i is not a gluon, the soft limit is simply not relevant 
(the matrix element is non singular), and the S functions can assume arbitrary values. 
This implies that, in such a case, wc may replace 2Ei/^ with 1 in the definition of 
dij. Along the same lines, one may always replace (3i and/or (3j with 1 in dij when 
particle i and/or j is massive; in doing so, dij will not be proportional any longer to 
ki-kj when as — bs, but this is perfectly acceptable. We have tested these options, 
and will comment on this point in sect. 7. 
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The S functions defined in eq. (5.8) are suitable for numerical implementation, 
since they are well behaved in the whole {n + l)-body phase space, and in particular 
in the soft and coUinear limits. In fact, after trivial algebra (also using eq. (5.14)) 
one gets 

Vd,,^l+ T ^+ E ^hkiizki). (5.15) 

Ij^j k^{i,j} or 

The Ei dependence can be analytically removed in the second term on the r.h.s of 
eq. (5.15), and therefore Vdij can be computed numerically both in the soft and 
coUinear limits, which are needed for the evaluations of the soft, coUinear, and soft- 
coUinear counterevents. 



5.3 {n + l)-body matrix elements 

In the previous section we have seen how to compute the S functions in eq. (4.29) in 
the various singular limits. The other major ingredient for the evaluation of eq. (4.29) 
is the damped matrix element 

M\]^'''^ = il-m,)C!M^-'^'^'K (5.16) 

We shall now give explicit, finite, expressions for M.\"'^^''^^ in the soft and coUinear 
limits. Using the notation introduced in eq. (5.7), we have 



where a is the relabeling of the indices introduced in eq. (5.7). Note that the eikonals 
on the r.h.s. of eq. (5.17) are finite at = 0, since the factor in the numerator is 
cancelled by that resulting from Ef in the denominator (see eq. (4.26)). Furthermore, 
in the case nij = 0, the coUinear limit of eq. (5.17) is also finite, since the factor 
1 — yij in the numerator cancels an identical factor in the denominator, when p = j 
OT I = j (see eq. (4.27)); this coUinear limit can thus be used to construct the soft- 
coUinear counterterm. The cancellations of the ^| and 1 — yij factors are carried 
out analytically, and eq. (5.17) is therefore well behaved numerically in the whole 
{n + l)-body phase space. 
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In the case = rrij — 0, we also need to compute the coUinear hmit. We have 
(see ref. [12]) 



lim A.;r''°' {r: {*}) - «i^^^a,j%<)>^""°' (^»-^ {*}) 



+ 9l '\^'^' Qx,i;,.fa)A^'"'°' (r; {*'^'}) , (5^8) 

where Zji has been defined in eq. (5.9), and the functions Qab*{z) and reduced matrix 
elements are given in app. D. The second term on the r.h.s. of eq. (5.18) 

vanishes when the integration over the azimuthal angle (pi is carried out. However, 
it is in general non zero pointwise in the phase space and, as eq. (5.18) shows, it 
actually has the same behaviour for yij 1 as the first term; it must therefore be 
included when constructing a local collinear counterterm. 
Using eqs. (4.26) and (4.27), eq. (5.18) becomes: 



to^Allr"> ('-; {*}) = |f7'f'x,.,(%.)Al'-°' {r'^-\ {*}) 

+ Ji«W.,fe.)AT'-°'(r;{t'^"}). (5.19) 

where, in analogy to eq. (4.44), we have defined 

Q,,.{z) = {l-z)Qab*{z). (5.20) 

It is apparent that the r.h.s. of eq. (5.19) can be computed numerically also in the 
soft limit, Ei (which is equivalent to Zji — > 1), owing to the regularizing factors 
(1 — ^) of eqs. (4.44) and (5.20). In the cases (Xi,Xj) = {g,g) and {Xi,Xj) = {q,g) 
(where g is a massless quark) the Altarelli-Parisi kernels have a soft singularity for 
Zji — > 0. However, this singularity is eventually damped, in the computation of the 
cross section, by the factor hij{zij) in eq. (5.8), and thus need not be subtracted. 



6. Optimization 

What we have done so far allows one to implement the FKS subtraction method in a 

computer code. We shall now tackle the problem of how to make this straightforward 
implementation more efficient. As discussed in sect. 3.2, the number of contributions 
to the sum on the r.h.s. of eq. (4.28) scales approximately as the square of the number 
of light partons in the final state. We shall show in the following that within the 
FKS formalism it is trivial to reduce this scaling property from a quadratic form 
to a linear function at most, and actually to a constant if one only increases the 
number of final-state gluons. We shall also discuss in this section the integration of 
the n-body contribution, and the implementation of multi-channel integration. 



24 



6.1 Reduction of the number of independent contributions 

Since the number of contributions to eq. (4.28) coincides with the number of elements 

in PfkS) we start by observing that the definition given in eq. (3.37) is snfRcient 
to include all singularities to be eventually subtracted, but is redundant. As an 
example, consider the case of Xj = g, X,- = g, with q any massless quark. According 
to eq. (3.37), both (i,j) and (j, i) belong to V^ks- On the other hand, the pairs 
(i,j) and (j, i) are responsible for the same divergence of the matrix element, the 
coUinear one at ki \\ kj. This divergence is subtracted twice in eq. (4.28), in the terms 
da^j'^^^ and daj"'^^\ but is not double counted thanks to the presence of hij{zij) and 
of hji{zji) in the relevant S functions (see eq. (5.14)). Furthermore, as discussed in 
sect. 4.2, when the FKS parton is massive, all of the contributions da^J^^^ are finite 
without subtraction; it is therefore not particularly advantageous to have a massive 
FKS parton. Several other examples can be given, the upshot of which is that it is 
more convenient to define Pfks as follows, rather than as in eq. (3.37): 



Equation (6.1) differs from eq. (3.37) in the upper limit for the range of i, which 
implies that FKS partons are now restricted to be massless (this, in turn, has the 
consequence that the the condition kj — > need not be used any longer) ; and because 
of the presence of non-redundancy conditions, which we specify as follows: 



These conditions fully remove the symmetry between and (j, i) (for i and j both 
in the final state), except in the case in which both of these partons are gluons; for all 
other flavour combinations, either of these pairs does not belong to T'fks according 
to the new deflnition. 

The non-redundancy conditions are easily understood. Equation (6.3) implies 
that for those (unordered) pairs of particles which do not induce soft singularities 
(such as the qq pair of the example above), of the two possible ordered pairs only one 
is included in Pfks- On the other hand, since only the soft singularities associated 
with FKS partons are subtracted (see eq. (4.29)), it is useless to give to particles that 
are not gluons the role of FKS partons if their sisters are gluons. Therefore, the only 
case in which one needs a symmetric role for the FKS parton and its sister is that in 
which both are gluons. It should be clear, however, that the a-symmetrization carried 
out in eqs. (6.1)-(6.3) need be accompanied by an analogous a-symmetrization in the 



7'FKs(r) = |3<^<nf^+2,n/<J<nf^+n^^ + 2,i7^j, 

_y^(n+l,0)^^ynf ) ^ OO if k^ ^ OT k \\ kj, 




(6.1) 



li = g, Ij 7^ G Pfks 

¥"9, ¥"9, ihj) e ^FKS 



(j,O^Pfks if 3<j, (6.2) 
(J,^)^PFKS if 3<i<i. (6.3) 
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definition of the S functions. This can be achieved by replacing eqs. (4.21) and (4.22) 
with the following definitions: 

h,,{z) = 1 a li^g or j^l or j = 2, (6.4) 

hij{z) = h{z) if Ti = g and Ij = g . (6.5) 

It is immediate to see that all formulae given in sects. 3, 4, and 5 still hold, pro- 
vided that one understands eqs. (6.1)-(6.5) where relevant, rather than eqs. (3.37), 
(4.21), and (4.22). In the rest of this paper, and in the implementation of the FKS 
subtraction in MadFKS, we always use eqs. (6.1)-(6.5). 

At this point, the sum on the r.h.s. of eq. (4.28) still scales quadratically with 
the number of light partons; however, that sum obviously contains several identical 
contributions. Suppose for example that one has m final-state gluons ii . . . im- There 
will be m{m — 1) gluon pairs in Pfks, which will all contribute to da^"'^^\ However, 
it is clear that da^"^^^^ — da[^^^^^ , because of the symmetries of the matrix elements, S 
functions, phase space, and (IR-safe) observables. Likewise, if j is not a gluon, and is 
such that G Pfks, then {ia,j) G Vfks, for 2 < a < m, and da^^^^^^ = d(xl^^^\ 

In general, the determination of which elements of Vfks give identical contributions 
depends on the process considered, and on the underlying theory (e.g., QED only 
or full EW interactions), and can be implemented with relative ease in a computer 
code. It is clear, however, that if several identical particles are present in the final 
state, one and only one of them (which one in particular, is irrelevant) will be a 
member of the FKS pairs. We shall denote by 

Vfks ^ ^fks , (6.6) 

the subset of T^fks whose elements give non-identical contributions to the sum in 
eq. (4.28). We shall then compute 

da("+^)(r)= J2 ^t"\r)da^^'\r)^ (6.7) 

with ^l^^^"* an integer symmetry factor, equal to the number of identical contributions 
(for a given (i, j')) to the sum in eq. (4.28). The number of elements in 7^fks('") and 
^FKs(^") will be denoted by A^fks(^") and A^fks(^") respectively; the value of A^fks for 
several production processes will be reported in sect. 7. 

6.2 n-body matrix elements 

The integration of the n-body contribution, eq. (4.3), is straightforward, and is usu- 
ally carried out independently of that of the {n -\- l)-body contribution. In this 
section, we give the formulae that allow one to integrate these two contributions at 
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the same time. For instance, this technique is the default in MC@NLO. We start by 
observing that 

if e TZn =^ 3r/3 e T^n+i such that = with Xi — g. (6.8) 

This states the fact that, given a set of particles that corresponds to a physical 2 ^ n 
process, by adding one gluon in the final state one gets a physically meaningful 
2 — > n + 1 process. It therefore follows that 

^ d(T^'^\ra) = J] d(T^'^\r\) for a given i with 2i = c/. (6.9) 

Note that, for an arbitrary 2 — > n + 1 process rp, the corresponding 2 ^ n process 
does not necessarily exist (e.g. if Ti = q). In such a case, we have (i(7*^"^(r^) = 0. 
On the other hand, if one fixes i, if exists then it is unique, and hence eq. (6.9) 
holds. 

Let us now introduce the shorthand notation 

= (6.10) 

The precise form of B can be read from eqs. (4.4), (4.5), (4.12), and (4.14), but is 
not relevant in what follows. Fixing i and assuming is a physical 2 — > n process, 
we have, using eq. (4.31), 

da(-) ^B-^ Q{Ut - ^i)d^idyi^d^d4^i{^i = 0) . (6.11) 

Scut "5 

Using eq. (4.18) we can also write 

= = 0) 7 ®(^-* - ^i)d^idyijd^id(pm = 0) • (6.12) 

This equation holds for any given i with Xi — g. It is convenient to give a prescription 
for fixing i given the process. In order to do that, we observe that there is only one 
gluon among the first elements of the pairs that belong to T^fksj and this gluon is 
therefore a natural candidate for fixing i. We rewrite eq. (6.12) as follows: 

da^^\r)^ E ^oiM^'ir^. ren^^^, (6.13) 

{i,j)€VFKS 

daif = BS,,{^i = 0) ^ e{Ut - ii)diidyi,d^d4^i{ii = 0) , (6.14) 

where on the l.h.s. of eq. (6.13) the use of a 2 — > n + 1 process as an argument for 
dcr^") is justified by the fact that there is a one-to-one correspondence between r and 
if one considers 'Pfks- 
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After these manipulations, the physical cross section of eq. (4.1) is written as 
follows: 

dap^P^{Ki, K2) = 

X (d(7("+^)(r; xi/fi, X2K2) + da("+^)(r; xi^i, X2K2) + da^''\r; xiKi, X2K2)) , 

(6.15) 

with (ia-("+^) given in eq. (4.40), and 

da(-+^)(r) + daH(r) = ^ + V.^^V^)) • (6-16) 

By inspection of eqs. (6.14) and (4.38), we see that the second term on the r.h.s. of 
eq. (6.16) has the very same structure as the soft counterevent of the first term. When 
cast in this form, it thus becomes obvious that the (n + l)- and n-body contributions 
to an NLO cross section can be integrated together. 

It is also easy to show that the degenerate (n + l)-body contribution da^^'^^^ can 
be integrated together with the other two contributions. This is the default strategy 
in MC@NLO. We refrain from reporting the relevant derivation here, and postpone 
it to a forthcoming paper, where the implementation of the FKS subtraction will be 
extended to the case of hadronic collisions. 

We conclude this section by observing that the logic behind eqs. (6.15) and (6.16) 
is that of getting the underlying Born process given a real-emission process. This 
structure matches quite smoothly that of MadGraph, which is the reason why we 
have chosen it in the first version of MadFKS. On the other hand, one may also 
consider the inverse logic, namely that of getting the real-emission process(es) given 
a Born process, which is especially convenient in the context of the matching of NLO 
computations with parton showers. We shall describe this option in app. E. 

6.3 Multi-channel integration 

In the previous sections we have shown how to define a number A'^fks of finite and 
independent short-distance contributions, whose sum is the observable cross section. 
For large final-state multiplicities, the integration of each of these contributions may 

prove difiicTiIt from the numerical point of view. 

Let us start by summarizing the strategy adopted in a tree-level computation by 
MadGraph/MadEvent [33]. If the tree-level matrix elements are those relevant to a 
2 — > n -|- 1 process, one can use eq. (4.1), set da^"'~^^^ = and da^'^^ = there, and 
define 

da^^+^\r) = M^''+''^\r) ... . , (6.17) 
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where J'^l '+i denotes that one reconstructs exactly n^^'* + 1 jets in the final state, 
which is sufficient to cut ofi^ all phase-space singularities of Let us now 

denote by 

da{r), a = l,...Nd{r) (6.18) 

the subset of Feynman diagrams that contribute to and that cannot be 

obtained from each other with a permutation of four-momenta. In other words, 

^(n+1,0) j-^-j _ ^da{r) -|- permutations of momenta^ . (6.19) 

a=l 

MadGraph then splits the integral of ^(7*^"+^^ in cq. (6.17) into Nd{r) independent 
contributions^^, called integration channels, and defined as follows: 

= Yl d(^a^^\r) , (6.20) 

a=l 

where 

d<^'r'\r) = ■^^""^•°nr)^#n+i, (6.21) 

D{r) = J2 \Mr)f ■ (6.22) 

/3=1 

The formal similarity between what done here and what is done in the context of the 
FKS subtraction formalism suggests to replace the damped (n + 1) matrix elements 
of eq. (5.16) with 



A^Sr'°Hr) = (l-y,,)e,^M("+^'°)(r) 



j(n+l,0) 
ij 

|2 



Mt,X'''\r) = (1 - mMM^''^'''\r)\-^^ , (6.23) 

and apply the subtraction procedure to this quantity. Unfortunately, this is not going 
to work. In the soft limit ^ 0. all diagrams that would induce a divergence in the 
square of the amplitude in eq. (6.19) vanish if squared individually, as in eq. (6.22), 
owing to the presence of a self-eikonal of a massless particle (an exception is therefore 
that of diagrams in which a gluon is emitted by a final-state massive particle). It 
follows that the cross sections da^^^^ that are equal to zero are actually those that 
one would hke to have the largest contributions to the sum in eq. (6.20). 

Prom the physics viewpoint, this is easy to understand. The squares of individual 
Feynman diagrams work well in tree-level computations because these computations 



^^We point out that MadGraph does not include diagrams that contain four-point vertices in the 
weights used for multichanneling. This fact is irrelevant for the sake of the present discussion. 
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are basically classical physics. On the other hand, soft singularities are inherently a 
quantum effect, since they are due to the interference of diagrams. Thus, multichan- 
nel sampling based on squares of diagrams cannot possibly describe all complications 
due to interference. 

We can bypass this problem in the following way. We first introduce a monoton- 
ically decreasing, smooth function f{x) with the following properties: 



< fix) < 1 
lim f{x) — 1 . 



for < a; < 1 , 



(6.24) 
(6.25) 



,(n+l,0). 



We then modify eq. (4.24) to read 

5,,(r)A^("+^'°)(r) 



r) 



{i,j)€PFKS 

+ E 

(«j)ePFKS 



1-/ 



(6.26) 



The factor in square brackets in the second term on the r.h.s. of eq. (6.26) damps the 
only soft and coUinear singularities that survive the damping of Sij. That second term 
is therefore finite over the whole (n + l)-body phase space, and is thus effectively 
a tree-level computation that can be dealt with in the same way as in standard 
MadGraph. 

The effect of the function / in the first term on the r.h.s. of eq. (6.26) is that 
of suppressing hard and large-angle emissions of the system. We can therefore 
assume that the typical kinematic configuration will always be well approximated by 
the degenerate configuration obtained by letting particles i and j become coUinear, 
or when particle i becomes soft. As we have seen in sect. 5.1, these two limits result 
in the same four-momentum configurations, and we can choose one or the other 
depending on which of the two will result in a physical 2 ^ n process as far as 
particle identities are concerned. In practice, if Xi = g we shall use the soft limit, 
while if Xi ^ g the coUinear limit will be used. To avoid unnecessarily complicating 
the notation, we shall denote both cases by in this discussion. Starting from 
eq. (6.26), we get 



(ri+l)/ \ r (n+l), 



(6.27) 



/3=1 



with 



Ml 

D{r) 



E 



1-/ 



5i,(r) lA<("+^'°)(r) 



(S) 

Jf{r) 



(6.28) 
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We stress again that, in spite of the presence of J'^l rather than that of J"^^ the 
cross section in eq. (6.28) is finite, and needs no subtraction. We have also defined 

0-^^^ = ^/ (^) "-S'^) . (6.29) 

where dcr^J^^' has been given in eq. (4.29), and we appUed in eq. (6.27) the same 
symmetry factor (^ij'^^^ as in eq. (6.7) for obvious reasons. As discussed above, the 
kinematic configuration used to evaluate the prefactor \da\'^ /D in eq. (6.29) is that 
of either the soft or coUinear counterevents, which coincide. Therefore, this prefactor 
will have the same value in the case of the event and of all the counterevents. On 
the other hand, f{2ki-kj/s) will in general be different when computed with event 
and counterevent kinematics, but the smoothness of this function guarantees that 
the cancellation of phase-space singularities will still take place. 

Note that the formulae above can be applied without modifications even if one 
sets f{x) = 1. This means that n-body Fcynman diagrams are used for multichannel 
sampling in the whole {n + l)-body phase space. This may imply a loss of efficiency 
when hard-emission regions are integrated over. However, as we shall see in sect. 7, 
this is not the case for any of the processes we have considered in this paper. 

From cq. (6.27), we can read the total number of integration channels of our 
implementation of an NLO cross section, for a given process r: 

A^ch(r) = Na{r) + J] iVrf(r^) . (6.30) 

{«J)€Pfks 

Note that the second term on the r.h.s. of this equation is roughly proportional to 
-^FKS (it would be exactly so if all Nd{r^) were equal). So while A'"fks is directly re- 
lated to the performance of the FKS formalism in keeping the number of independent 
subtraction terms as small as possible, the ratio iVch/A^FKS can be used to measure 
the increase of the number of integration channels due to the inherent complexity 
of the underlying Feynman diagram structure. Note that, by setting f{x) = 1, 
eq. (6.30) becomes 

iV,h(r) = Yl ^'^(^'') • (6-31) 

7. Results 

In this section we present MadFKS results for several production processes in e+e" or 
coUisions. These correspond to the full NLO cross section, eqs. (4.1) or (6.15), 
where we set 

V^Fiii = (7.1) 

in eq. (4.14). This is equivalent to assuming that the one-loop contribution is a 
pure-pole term in the CDR scheme (see app. B), which implies that the results 
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given in this section are non physical. Indeed, our aim is not that of carrying out 
a phenomenological study, but to document the convergence properties of the FKS 
subtraction in MadFKS (i.e., the statistical errors obtained with a given number of 
integration points), and to make sure that our implementation is correct. 

We point out that the optimizations implemented so far are restricted to the FKS 
formalism; from the point of view of MadGraph, the current version of MadFKS must 
be considered as a benchmark upon which we shall build a fully efficient program. In 
particular, there are two main lines of development for the code. Firstly, the results 
we present here have been obtained by running all integration channels with the 
same number of points. This procedure will be improved following the same strategy 
as is used in MadGraph, namely that of performing preliminary low-statistics runs 
to determine which channels give significant contributions, in order to further the 
runs with larger statistics only for those. Secondly, in the current version of MadFKS 
all sums over helicities are performed exactly. A much more efficient procedure is 
that of performing such a sum using Monte Carlo techniques, as is typically done in 
the context of tree-level computations. This is also possible at the NLO, provided 
one is able to write the subtraction formalism for any given helicity configuration. 
This is straightforward in the context of the FKS subtraction method, since indeed 
all necessary ingredients arc already available in the literature [12, 34]; we shall 
further comment on this point, and give all the relevant details, in app. F. We have 
refrained from implementing these improvements (and several other minor ones) in 
the first version of MadFKS, since this is the best way to gauge in a clear manner the 
capabilities of the code; for example, by treating all channels on the same footing, 
one is able to check that they are all computed correctly, regardless of whether some 
are numerically relevant or not. 

As explained in sect. 6, to obtain a cross section we find it convenient to sum 
over (n+ l)-body processes, according to eq. (6.15). In this section, therefore, when 
considering a given partonic process, it will always be understood as an (n + 1)- 
body one. The corresponding n-body contributions will be added according to the 
procedure described in sect. 6.2. 

We have set up two different types of numerical tests to check the correctness of 
the implementation. The first type of tests make sure that the local behaviour of our 
main ingredients in the computation of the real-emission contribution, namely the 
matrix elements and the S functions, behave as we expect them to do in the coUinear 
and soft limits. For all real-emission processes r, and for all pairs G Vfks, we 
check numerically that eqs. (5.17), (5.18), and (4.16)-(4.20) arc satisfied. This is done 
automatically, and if one of these tests fails, the programme proceeds no further. 

The second type of tests rely on the fact that several free parameters are intro- 
duced in the FKS short-distance cross sections; the final results must be independent 
of these parameters, but partial results must not. To give one example, let us con- 
sider the case of ^cut, introduced in sect. 4. This parameter enters analytically the 
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n-body cross sections (see e.g. eq. (4.6) and the results for the integrated eikonals in 
app. A), and numerically the (n + l)-body cross section through the definition of the 
soft counterterm, according to eq. (4.33). Therefore, the quantities (icr^""*"^^ and rfcr^^^ 
will separately depend on ^cuti but their sum will not. The free parameters whose 
variations we have considered here are ^cut, ^Oi cis, bs- We have also investigated 
if the final results are independent of whether one uses the factor lE^I ^ in the 
definition of dij (this option is denoted by useenergy=.true.), or one uses 1 instead 
(this option is denoted by useenergy= . false . ); see sect. 5.2 for more details. Each 
of these parameter choices corresponds to a set of (parallel) runs. Therefore, to be 
definite, we have used the [n + l)-body process 

e+(l) e-(2) ^Z^ u{^) fZ(4) ^(5) ^(6) ^(7) (7.2) 

as a benchmark; this is a contribution to the four-jet cross section at the NLO. We 
have run at ^/s — 100 GeV, and set the factorization, renormalization and Ellis- 
Sexton scales equal to the Z mass, /J'j? = fJ^R = = Tn\. This process is indeed 
simple enough to run in a short amount of time, but has the required complexity to 
test the implementation in all its details. The condition in eq. (3.36), i.e. the hard 
cuts, has been imposed with the KTCLUS routine that implements the jet-finding 
algorithm of ref. [35], with Fcut = (10 GeV)^. The jet four-momenta are defined as 
the sum of the four-momenta of the partons in the jets; the lowest jet energy is 
therefore about 10 GeV. The cross sections and integration errors (as returned by 
Vegas [36]) for the process in eq. (7.2) are given in table 1. There are A'"fks = 3 
independent FKS pairs contributing to this process (these can be chosen to be (6, 3), 
(6,4), and (6,7) according to the labels of cq. (7.2)), and five integration channels 
for each of them. We used 50000 integration points and 10 iterations per channel; 
only about 15% of these points result in kinematic configurations that pass the jet- 
finding conditions; at the level of phase space, the points are generated flat, i.e. no 
information on the hard cuts is available when one generates the kinematics. The 
results of table 1 are clearly independent^^ of the free parameters, typically within 
one standard deviation. As an even tighter consistency check, we have picked the 
largest and the smallest values among those in the table, and re-evaluated them with 
higher statistics. For a six-fold increase of the number of iterations, keeping the 
points per iteration to the same value as before, the corresponding cross sections 
change as follow: 

3.5633 ±0.0154 ^ 3.6086 ± 0.0051 , , 

7.3 

3.6350 ±0.0151 ^ 3.6007 ±0.0053. 
It is interesting to observe that, apart from the convergence of the mean values 
towards a common value, the errors in eq. (7.3) scale approximately as we would 

-'^^Sincc in this specific example we do not sum ovqi flavours, we have to set Nf = in eqs. (4.7) 
and (4.9); this corresponds to not allowing g qq splittings at the Born level, which matches the 
fact that we do not include any (n -|- l)-body matrix elements other than that in eq. (7.2). 
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as — bs ^cut — ^max ^cut — 0.3 ^c«t — 0.1 ^c«t — 0.01 



useenergy= . true . 
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.5966 ± 0, 
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3, 


.5938 ± 


0, 
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3, 
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lb 
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3, 


.6053 ± 
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3, 


.5998 ± 0, 


,0181 


3, 


.5988 lb 


0, 


,0122 


3, 


.6088 


lb 


0, 


.0165 


0.06 


1, 


.5 


3, 


.6144 ± 





.0161 


3, 


.5986 ± 0, 


,0140 


3, 


.5847 lb 


0, 


,0119 


3, 


.5884 


± 





.0126 




2, 


.0 


3, 


.5990 ± 





.0166 


3, 


.6016 ± 0, 


,0158 


3, 


.6014 ± 


0. 


,0147 


3, 


.6191 


zb 





.0133 
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3, 
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3, 
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0, 


,0165 


3, 


.5986 
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.0142 




2, 


.0 


3, 


.5921 ± 





.0125 


3, 


,6260 ± 0, 


,0211 


3, 


.6034 lb 


0, 


,0134 


3, 


.6007 


zb 
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0.6 
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.5944 ± 0. 
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3, 
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0. 


,0123 


3, 
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3, 
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0. 
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0. 
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.0 
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3, 


.5942 ± 0. 
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3, 
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0. 
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3, 
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0.06 


1, 
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3, 
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3, 


.5977 ± 0. 
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Table 1: Cross section (in pb) and Monte Carlo integration errors for the (n + l)-body 
process e+e~ — Z — uuggg. See the text for details. 



expect them to do if they would follow the Gaussian law typical of integrals of 
ordinary functions (while a subtracted cross section is actually a distribution). This 
gives us confidence on the fact that Vegas estimates correctly the integration errors. 
We have also checked that the cross section is independent of the choice of the 
function h (see eqs. (6.5) and (4.23)); since the numerical effects are even smaller 
than those reported in table 1, we refrain from presenting the corresponding results 
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here. We have performed the runs with 



(1 - x)'^"- 

f{x)= a = l, c=l. (7.4) 

It is remarkable that, with this choice, the contribution of da^^^^^ (see eq. (6.28)) is 
smaller than the statistical errors reported in table 1, and is therefore not included 
in the results of the table. 

As a final comment on the process in eq. (7.2), we point out a fairly pleasant 
feature that could have been anticipated. Namely, the results presented in table 1 
have integration errors that are by far and large independent of the choices of the 
free parameters, thus including the case of "maximal" subtraction, ^cut = Cmax and 
= 2, when for each event one always subtracts the three counterevents. This is 
quite an useful property, since it implies that it will not be necessary to scan the 
space of the free parameters in order to find those values that maximize the numerical 
stability, all choices being almost equivalent. At the core of this property stands 
the fact that, for any given FKS pair (i.e., a finite and independent contribution), 
when scanning the (n + l)-body phase space one finds at most one soft and one 
coUinear singularities which, as discussed in sect. 5.1, are associated with a single 
kinematic configuration. Therefore, the behaviour of the subtracted cross section 
will not change much by varying ^cut and 5. It is instructive to compare this with 
the case of dipole subtraction. There, by changing the parameter that is typically 
called a (the analogue of our ^cut and 5), one can change drastically the number of 
subtraction terms included in the computation, and the tuning of a may become 
necessary. The capability of the FKS subtraction method to give numerically stable 
results also for large values of ^cut and 5 is especially promising in view of the fact 
that such values are usually the optimal choices when matching NLO computations 
with parton showers in MC@NLO or POWHEG. 



The independence of the results of the parameters ^cut and 5o allows one to 
check the correctness of the subtraction terms for the real-emission matrix elements, 
and of the logarithmic terms in the n-body contributions. In order to also check the 
non-logarithmic terms that enter the n-body contributions, we have computed the 
physical NLO cross sections for the processes e+e~ ^ Z ^ 2 jets and H ^ 2 jets, 
where if is a SM Higgs with rriH = 120 GcV, using the readily available relevant 
matrix elements for the virtual corrections. We have found complete agreement with 
the known results. We point out that, by doing so, we have tested all analytical 
formulae presented in this paper, except for eqs. (4.41) and (4.42), which only enter 
hadronic cross sections, and for the integrals of the non-massless eikonals, eqs. (A. 8), 
(A. 10), and (A. 12). The latter integrals, however, have been checked by computing 
them numerically with very high precision. 
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(n + l)-body process 



cross section (pb) 



Nfks 



e~^e~ — > Z ^ uugg 
e'^e" — > Z — > uuggg 
e^e~ — > Z — > uugggg 



1*/Z ^ jjjj 
1*/Z ^ jjjjj 
e+e" ^ Y/Z 333333 



e+e 



e^e — > Z — > ttgg 
e^e~ ^ Z — > ttggg 
e+e" ^ Z ^ tigggg 



e+e 



e^e 



Z ttbbg 
Z — > ttbbgg 



e+e ^ Z ^ ^i^tS'S'S' 
e+e- ^ Z ^ ggggg 



H gggg 
H ggggg 
H gggggg 



0.4144 ± 0.0006 
0.3601 ± 0.0014 
0.8869 ± 0.0054 



0.1801 ± 0.0002 
0.1529 ± 0.0004 
0.3954 ± 0.0015 



0.1219 ± 0.0003 
0.1521 ± 0.0013 
0.1108 ± 0.0031 



0.1972 ± 0.0024 
0.2157 ± 0.0029 



0.3712 ± 0.0037 
0.1584 ± 0.0019 



0.1404 ± 0.0005 

0.2575 ± 0.0018 
0.1186 ± 0.0008 



0.15% 
0.38% 
0.61% 



0.12% 
0.26% 
0.38% 





3 


))xlO^ 


3 


))xlO-^ 


3 


))xl03 


14 


))xl02 


30 


))xlO° 


55 



0.24% 
0.83% 
2.76% 



1.23% 
1.34% 



1.00% 
1.23% 



0.34% 

0.69% 
0.70% 



xlO-i 3 
xlO-2 3 
xlO-3 3 



xlO-^ 
xlO-^ 



xlO-s 
xlO-i 



4 
5 



3 
2 



xlO-^ 1 

xlO-s 1 
xlO-9 1 



Table 2: Cross sections and Monte Carlo integration errors (in absolute value and as a 
fraction relative to the cross section) for various processes. 



We have considered a variety of other processes, in order to prove the flexibility of 
the code, and to check the convergence of the numerical integration; the processes and 
the corresponding results are reported in table 2. We stress that, in order to obtain 

the results displayed in that table, we just had to write the relevant process in an 
input card for MadFKS; this is the same procedure as that done in MadGraph. All cross 
sections appearing in table 2 have been computed by setting ^cut = 0.1, 6o = 0.6, 
as = 1.5, bg = 1.5, and useenergy=.true. . We have run 10 Vegas iterations; 
for each of them, phase-space points have been generated in order to obtain about 
10000 kinematic configurations per iteration that pass the jet-finding cuts, which 
corresponds to the statistics used to generate the process in eq. (7.2). For final 
states with only massless partons, we set y/s — 100 GeV; we have used the notation 
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j to indicate the sum over partons (we have considered four flavours here); thus, 
e.g. e+e^ jjjjjj iSj up to the missing flnite one-loop contribution, the physical 
five-jet cross section. For top-quark production processes, we set ^/s = 500 GeV. 
The b quark is taken to be massive {mi, = 4.7 GeV) in processes with ttbb + X 
final states, and it does not enter the jet-finding algorithm. Finally, for stop and 
gluino production we have used m^^ = 400 GeV, rrig — 400 GeV, and ^/s — 1 TeV; 
we assume an effective Zgg axial vertex, and set the corresponding coupling equal 
to one. We have also studied the production of fully gluonic final states arising 
from the decay of an (off-shell) SM Higgs with rnn = 120 GeV, that contribute to 
three-, four-, and five-jet observables. In order to keep this process on the same 
footing as the others that feature only massless final-state particles, we have set 
y/s = 100 GeV, but used a //.'*' initial state in view of the smallness of the electron 
Yukawa coupling. As one can see from table 2, fuUy-gluonic Higgs decays have a 
single non-trivial contribution from FKS pairs. In the cases of processes involving 
massive particles we also tested the option of replacing (3^ and/or /3; with 1 in the 
expressions of the relevant dki (see sect. 5.2). We did not observe any significant 
differences with respect to the results obtained with the default choice. 

As one can see from the table, the numerical errors are fairly small even with 
the limited statistics we used. The largest of them (but still rather modest) are 
associated with the processes with a ti pair in the final state, plus either four gluons 
or a bb pair and extra gluons. It is clear that these processes feature several mass 
scales quite different from each other - the cm. energy, the top mass, the b mass, 
and the minimum jet energy (~ 10 GeV) - and therefore one expects the coefficients 
of the perturbative series to be plagued by several large logarithms, that inherently 
increase the complexity of the numerical calculation. We point out that, since we 
treat the b quark as massive, the bb pair is not an FKS pair, according to eq. (6.1). On 
the other hand, given the smallness of the b mass, it may be beneficial to treat such 
a pair as an FKS pair, which would also possibly imply the definition of a "quasi- 
coUinear" counterterm, whose analytical form has to tend to that in eq. (5.18) in the 
limit nib 0. We leave the implementation of this option to a future work. 

It is important to note that the growth of A'"fks with the final-state multiplicity 
is always rather modest. From table 2, one can actually see that A'"fks is a constant 
if the number of gluons is increased (A'"fks for ttbbgg is one unit larger than in the 
case of ttbbg, since only in the former process one can form an FKS pair with two 
gluons). The increase of A^fks in the e~^e~ — > (n + l)j processes is mainly due to 
the corresponding increase in partonic subprocesses: in the case of four, five, and six 
particles in the final state of the (n + l)-body processes, the number of contributing 
subprocesses is equal to 7, 7, and 17 respectively. Among individual contributions to 
the five-jet cross section, those with the largest A'"fks are the four-quark, different- 
fiavour processes such as e"^e~ — > uuddgg, which have seven possible FKS pairs (since 
one distinguishes between quarks and antiquarks when pairing them with gluons, in 
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Figure 1: Differential spectra for the first three partonic processes listed in table 2. 
The histograms for the latter two processes have been rescaled (by a factor of 9 and 250 
respectively) in order to fit into the layout. We present thrust, C parameter, and the 
energy and polar angle of the leading jet. 



order not to neglect possible charge asymmetries). 

As discussed in sect. 6.2, we integrate the n-body matrix elements at the same 
time as the (n + l)-body ones. On an event-by-event basis, we can therefore obtain 
both the NLO and the LO contributions. We have checked that the latter is, for 
all processes, fully consistent with the one predicted by standard MadGraph. If one 
switches off in MadGraph the optimizations relevant to the separate treatment of 
different integration channels, our LO computation has the same statistical accuracy 
as that in standard MadGraph. More importantly, if we only integrate the Born 
contributions to the processes listed in table 2 with the same number of points as 
that used for the NLO contributions (distributed equally among the possibly smaller 
number of integration channels), the resulting integration uncertainties are a relative 
factor 1.9 to 4.5 smaller than those relevant to the NLO results presented here. 
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E.„d (GeV) cose,^^^ 

Figure 2: Same as in fig. 1, for the three- (black sohd), four- (red dashed), and five-jet 
(blue dotted) cross sections (fourth to sixth processes in table 2). 



Exceptions are found for the (n + l)-body processes tibbg and tigggg, whose LO 
contributions have integration accuracies better than the NLO ones by a factor 7 
and 9 respectively. Overall, these figures give us another indication of the fact that 
the subtraction of singularities is achieved in very satisfactory manner from the 
numerical point of view. 

We now turn to presenting differential distributions for some of the processes 
listed in table 2. As clarified in sect. 5 (see in particular sect. 5.1), for a given process 
and choice of E Vfks, to each random number there correspond two kinematic 
configurations (for the event and for the counterevents plus n-body contributions), 
with associated weights. These configurations give the complete information on all 
the four-momenta of the (resolved) final-state particles, and thus one can plot as 
many (infrared-safe) observables as one wants in the course of a run. The plots we 
present here are obtained with the same statistics (in fact, in the same runs) as that 
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Figure 3: Differential spectra for tiie seventii to nintli partonic processes listed in table 2. 
The histograms for the latter two processes have been rescaled (by a factor of 6 and 50 
respectively) in order to fit into the layout. We present the energy and polar angle of the 
top quark and of the leading jet. 



used for the total cross section results given in table 2. We point out that we do 
not use any smoothing procedure for the spectra we show: we limit ourselves to 
filling the histograms putting the weights in the bins determined by the kinematic 
configurations given by the code. This implies, in particular, that large weights for 
the event and counterevents may fall into different bins. We choose a relatively fine 
binning, in order to expose in a clear manner whether this mis-binning (which is 
unavoidable in any subtraction-based computation) is a severe problem or not. 

In fig. 1 we show the distributions in the thrust, C-parameter, and energy and 
polar angle of the hardest jet, for the first three processes listed in table 2. In 
fig. 2, the same distributions are shown for three-, four-, and five-jet final states, and 
including the Z/7 interference effects - the plots in fig. 2 would therefore correspond 
to physical jet cross sections, had we included the proper one-loop contributions. 
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Finally, in fig. 3 we show the energies and polar angles of the top quark and of the 
hardest jet relevant to the e+e^ ^ ti + {n + l)g processes. All spectra are fairly 
smooth, except those relevant to the polar angles, where fluctuations are marginal 
for light-jet-only production processes, and more evident in the case of tt processes. 

8. Conclusions 

In this paper, we have presented a complete automation of the computation of any 
cross section at the next-to-leading order in QCD, with the exception of an infrared- 
and ultraviolet-finite term of one-loop origin. This is achieved by embedding into 
the MadGraph framework the universal subtraction formalism of Frixione, Kunszt, 
and Signer (FKS); the resulting computer code has been named MadFKS. The only 
intervention needed to run the code is the definition of the inputs, among which 
is the production process one wants to study. The process is any physical process 
that results from a theory implemented in MadGraph, thus including the Standard 
Model and any user-defined Beyond the Standard Model theory. The current version 
of the implementation allows one to compute cross sections with colourless particles 
in the initial state, but there are no difficulties in extending it to other types of 
colliding particles; we defer this to future work. It is also worth mentioning that the 
formalism presented here requires only trivial modifications to be extended to the 
subtraction of QED infrared singularities. With minimal effort, one could therefore 
fully automate the FKS subtraction for the complete electro-weak corrections; we 
point out, however, that in order to actually perform such computations it will 
probably be necessary to implement a complex-mass scheme into MadGraph. 

The key features of our implementation are a direct consequence of the FKS sub- 
traction formalism. Through a dynamic partition of the phase space, one effectively 
defines partonic processes which are either non singular, or have at most one soft 
and one collinear singularity, thus rendering the subtraction of these singularities a 
fairly easy, and numerically efficient, task. These partonic processes are separately 
finite and fully independent of each other, and therefore the method is naturally 
suited to parallel computations. For each phase-space point, there are only two inde- 
pendent kinematic configurations; one corresponds to a fully-resolved configuration 
(all partons are hard and well separated from each other), and the other one is the 
associated partly unresolved configuration (where either one parton is soft, or two 
partons are collinear, or both), whose weight is the sum of all subtraction terms, of 
all finite remainders of the subtraction procedure, and of the Born; this structure is 
exact, i.e. does not imply any approximation. Being able to limit the independent 
kinematic configurations to two is quite beneficial for the numerical stability of the 
results, and in particular for the smoothness of the predictions for differential spec- 
tra. The latter point is easy to understand, since in FKS one renders it small the 
probability that weights large in absolute value and opposite in sign will not fall in 
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the same histogram bin (thus, one reduces the mis-binning problem that affects all 
computations based on a subtraction method). 

The FKS subtraction method allows one to introduce several free parameters in 
the calculation; partial results (e.g., the (n+ l)-body and the n-body contributions) 
depend on the values chosen for these parameters, but the physical cross section does 
not. We have exploited this property to check the correctness of our implementation. 
In addition, we have shown that the integration errors are largely independent of 
the choices of the free parameters, which is a direct consequence of having only 
two independent kinematic configurations for each phase space point. This is an 
important fact, since it allows one to choose the free parameters without having to 
worry about the numerical stability of the results. 

We have presented predictions for several production processes in e~^e^ or 
collisions. The numerical integrations converge quickly, requiring relatively modest 
statistics to achieve accurate results. We have run in parallel, with each of the 
parallel jobs corresponding to one integration channel, and using the same statistics. 
This is certainly not the best way to perform a multi-channel computation. The 
most obvious improvement, which is straightforward to automate, is that of starting 
with a low statistics run for all channels, checking which contributions are largest 
in absolute value or/and have the largest fractional errors, and further the run with 
larger statistics only for those channels. We plan to release a public web interface to 
MadFKS as soon as the case of hadronic collisions will also be included, and some of 
the optimizations discussed in this paper will be implemented. 

The current version of MadGraph, upon which MadFKS is constructed, is based 
on the evaluation of Feynman diagrams. We point out, however, that MadFKS can be 
rendered fully independent of Feynman diagram information, thus implying the pos- 
sibility of using recursively-generated amplitudes. In fact, the only part in the present 
implementation which is strictly dependent on Feynman graphs is that relevant to 
multi-channel sampling. A different criterion for sampling, based e.g. on topological 
information or on the computation of the amplitudes, can be easily incorporated in 
MadFKS (where multi-channel sampling is the very last step of the computation) . 

If combined with one-loop computations, MadFKS will allow one to predict any 
NLO cross section, limited only by CPU time. The present bound is actually set by 
the capabihty of MadGraph to handle the large number of Feynman diagrams and 
colour configurations for very large final-state multiplicities; tree-level amplitude 
computations based on recursion relations will push this bound to larger multiplic- 
ities. The computer programme can now use any external routine that returns the 
finite one-loop contribution as defined in this paper, feeding it with four-momenta 
and information on the tree-level structure. It is clearly desirable to fully incorporate 
one-loop computations into the programme, which appears to be feasible given the 
recent progress made in the evaluation of virtual graphs. Either way, the implemen- 
tation presented here is the first step towards a full automation of NLO computations 
matched with partons showers, as done in MC@NLO or POWHEG. 
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A. Eikonal integrals 

We define the integrals of the eikonal factors as follows: 



with 

In eq. (A.l), Ski collects by definition all divergent terms (times an e-dependent 
pre-factor, to be given explicitly below), while Ski is finite. The normalization of 
eq. (A.l) is conventional, and is equal to that used in ref. [12], times a factor Stt^ 
(this factor has been compensated by a different normalization of the colour-linked 
Horn's M^^f^^ w.r.t. that of ref. [12]). The energy of parton i and the measure (in 
3 — 2e dimensions) over its angular variables dfli are defined in the cm. frame of 
the colliding partons (see eq. (4.25)). In ref. [12] only the massless case = and 
mi = had been considered. The cases rrik = 0, mi ^ 0, and m^ ^ with k = I, 
have been dealt with in ref. [31], but the analytical results for eq. (A.l) have not 
been published there. In this paper, we have computed the only missing ingredient, 
i.e. mfc 7^ and mi ^ 0. We summarize here all the relevant results. 

• mk — 0, I — k (massless self-eikonal) . 



4^ = 0, (A.3) 
4^ = 0. (A.4) 

These equations trivially follow from the mass-shell condition /c^ = in eq. (A. 2), and 
serve the sole purpose of excluding the contributions of M.^kk^^ with nj < k < ri~^^ + 2 
from eq. (4.12). 
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• mjfc — 0, mi — 0, I ^ k. 
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The definition of /3; has been given in eq. (5.12). 
• mk 0, I = k (massive self-eikonal) . 
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and defined 
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(A. 14) 



(A.15) 



A 



OikiEk — El , 



(A.16) 




(A. 17) 



It is straightforward to clieck tliat, if ki — k^, eqs. (A. 11) and (A. 12) coincide with 
eqs. (A. 9) and (A. 10) respectively. 

B. Finite one-loop contribution 

Assuming the vahdity of KLN theorem, the most general form of the divergent part 
of the ultraviolet-renormalized one-loop contribution defined in eq. (3.25) can 

be obtained by computing the divergent contributions resulting from real-emission 
diagrams, and changing the signs in front of the poles in e obtained in this way. These 
divergences may have coUinear and/or soft origin. The results for the former have 
been computed in ref. [12] (see sects. 4.3 and 4.4 there). The structure of the latter 
is more complicated, since it arises from the integrals of the eikonal factors, which 
in turn depend on the masses of the two particles entering the eikonal. A summary 
of the relevant results is given in sect. A; here, we shall need the divergent parts, to 
be found in eqs. (A.5), (A.7), (A.9), and (A.ll). Using eqs. (3.33) and (3.34), we 



obtain: 




(B.l) 



with 





k=ni 



l=k+l 





(B.2) 



45 



Equation (B.2) generalizes eq. (3.2) of ref. [12] to the case of arbitrary particle masses, 
and is consistent with the results given in ref. [37]. We stress again that the poles in 
eq. (B.2) have infrared origin, since the ultraviolet divergences are assumed to have 
been eliminated through renormalization. 

The Ellis-Sexton scale Q, originally introduced in ref. [32], is any scale that may 
be used in one-loop computations to express the arguments of all the logarithms 
appearing there as s^/Q^ rather than as Sy/ski, where and Ski are two invariants 
constructed with the four-momenta of the particles that enter the hard process (see 
e.g. sect. 6 of ref. [38]). Clearly, it is not mandatory to introduce the Ellis-Sexton 
scale in the computations of the virtual amplitudes, and fairly often one just uses the 
factorization or renormalization scale instead. On the other hand, if Q is kept distinct 
from both /ip and /ir, one sees that the NLO cross section is exactly independent 
of its choice; in other words, variations of Q do not lead to any estimate of missing 
higher orders in perturbation theory, but allow one to check the internal consistency 
of the implementation. Although indeed most of the available one-loop results do 
not use the Ellis-Sexton scale, it is always possible to introduce it; we show how this 
can be done in app. C. 

Once the virtual amplitude is computed, one can interfere it with the Born am- 
plitude and obtain A4'^"-'^\ Thus, eqs. (B.l) and (B.2) can be solved for the finite 
contribution Vp*"^^ needed in eq. (4.14). This operation can be unambiguously car- 
ried out only after specifying a scheme. Two typical situations may be considered: 
the Conventional Dimensional Regularization (CDR) scheme, and the Dimensional 
Reduction (DR) scheme. In the former, the quantities A^^"'°^ and M^i'^^ in eq. (B.2) 
are evaluated in 4 — 2e dimensions, while in the latter they are evaluated in 4 dimen- 
sions. The finite parts Vp'j^J defined in the two schemes can be easily related (see 
e.g. ref. [38]): 

VifJ (CDR) = V^'ii (DR) - ^ 7(Xfc) , (B.3) 

k=ni 

with 

7(g) = -C^, l{g) = -CA. (B.4) 

Note that the sum on the r.h.s. of eq. (B.3) is extended to light quarks and gluons 
only. Consistently with the fact that the quantities in eq. (B.3) are ultraviolet- 
finite, we understand that all scheme changes relevant to ultraviolet divergences 
have already been performed; this implies that, from the ultraviolet point of view, 
the term vJry^^(DR) is given in the CDR scheme. 

We stress again that V^^^ which is used in eq. (4.14) must be computed in CDR. 

C. How to set lip /iR 

The formulae presented in this paper assume /i — /ip — /ir, and that another 
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arbitrary mass scale, the Ellis-Sexton scale Q, is introduced in the computation 
of the one-loop corrections. In this appendix we explain how to relax the former 
condition, and how to insert the dependence on the Ellis-Sexton scale in an one- 
loop result that was computed without introducing it. We perform these tasks by 
imposing renormalization group invariance w.r.t. and ^r. We start by giving our 
convention for the scales that appear in the various formulae of this paper: 

If fip ^ fiR, all the formulae for the short-distance cross sections given 
in this paper must be computed with n = /ip, except for the argument of 
as, which must be set equal to fiR. 

The statement above implies in particular that the PDFs are computed at /i — /ip, 
as customary. The convention given above cannot possibly allow us to recover terms 
proportional to logarithms of ratios of scales, which we now proceed to determine. 
In order to do so, we rewrite eq. (4.1) or eq. (6.15) symbolically as follows 

dap,p, = /(^^) ^ /(^^) ^ f + + rfa^") + Clog 4) , (CI) 

with C the unknown quantity that we need to determine. In order to proceed, we 
have to specify the power of as that enters our cross section formulae; we denote by 
b an integer such that 

Note that in general we may have 6 = (e.g. in the case of a purely-EW process 
at the Born level). We can now impose the invariance of the physical cross section 
w.r.t. to the renormalization scale 

^ dap,p,^0, (C.3) 



where as usual this equation holds up to terms of O (ag"^^). Using the explicit forms 
of the short-distance cross sections that appear on the r.h.s. of eq. (C.l), and the 
convention stated above, we find that eq. (C.3) is equivalent to 



d 



frfa(^'") + Clog4') =0. (C.4) 



Therefore, since 



-^^^ - -Poasil^n) + 0{as) , Po ^2^ , (C.5) 

we obtain 

C = -27r/3o b (^^) c^^^^'") , (C.6) 
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where rfa*^^'") is the Born cross section, given in eq. (4.4). 

This would be the end of the story if all short-distance cross sections could be 
computed according to the convention given above. This may not be the case for 
the one-loop contribution. We can however always manage to rewrite the one-loop 
contribution using our convention, regardless of its original form. To do this, we im- 
pose the condition that the physical cross section be independent of the factorization 
scale: 

^ -dap,p,^0, (C.7) 



which is the analogue of eq. (C.3). After some algebra (where one also uses the 
evolution equations for the PDFs), eq. (C.7) is found to be equivalent to the condition 

^_^da^v,n) ^ _c , (C.8) 

with C given in eq. (C.6). In order to proceed, let us thus consider eq. (4.14), 
and denote by tk all possible quantities with mass-squared dimensions that can be 
constructed with four-momenta. Following our convention, we must write 

vi"yl)(/xL/^l,Q') = c.U^^l)WI^i^^%Q'), (c.9) 
v&lrV^, Q') = av {{tk}) log ^ + bv (g^ {t,}) , (c.io) 

where use has been made of the fact that, in an one-loop computation, there will be 
an explicit linear dependence on the logarithm of the renormalization scale; according 
to our convention, in the argument of this logarithm /ir must be replaced by fip, 
hence eq. (C.IO). On the other hand, the computation will in general contain terms 
of the kind log^t^/t^, with p = 1,2. The procedure of Ellis and Sexton implies that 
these must be replaced by (logt^/Q^ — logti/Q'^y, and by in eq. (C.IO) may thus 
include up to double logarithms of Q^; this is however unimportant in what follows. 
Equation (C.8) can now be seen as a relation to solve for ay] we obtain 

vfci(/^kQ^) = 27r/?o6 log^ + vi"yi^(g^g^) . (c.n) 

Note that the r.h.s. of this equation is indeed independent of /^r, since Ai^"''^^ is 
proportional to aK/x^j). Equations (C.9) and (C.ll) can be used to insert the renor- 
malization, factorization, and Ellis-Sexton scale dependences into a one-loop result 
originally computed with a single scale M. In other words, if we are given the 
quantity v{M'^) such that 

V^FmiM^ M\ M') = v{M^) , (C.12) 

then we can construct the finite contribution of one-loop origin which complies with 
our convention as follows: 

(/.L /^l, Q') = 2nPobM^-''^ log ^ + 4(/^l)^^ ■ (C.13) 
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We stress that the quantity v{Q'^) is proportional to q;|((5^), and therefore the whole 
r.h.s. of eq. (C.13) is proportional to ag{fij^). 

In summary, we shall compute the cross section with the most general assignment 
of scales IJ>f ^ I^r Q using 

dap.p, = + + - 27r/3o b ^^^^'"^ log ^) , 

(C.14) 

where the scales in the formulae for the short-distance cross sections that appear on 
the r.h.s. of this equation have to be chosen according to the convention given at 
the beginning of this appendix. When following this convention, the finite part of 
the one-loop contribution, eq. (4.14), can be derived using eq. (C.13) from a result 
t'(M^) obtained with a single scale M. 



D. Azimuthal terms in collinear limits 

When computing the collinear limit of the damped real matrix elements, we obtain 
eq. (5.18) (or the equivalent form eq. (5.19)). It is customary to write such a limit 
by only keeping the first term on the r.h.s. of eq. (5.18) but, as discussed in sect. 5.3, 
the second term is equally important if one wants to construct local countcrtcrms. 
Appendix B of ref. [12] explains in detail how to compute the reduced matrix clement 
and the kernels Qab*{.z). These kernels arc universal in the same sense as the 
Altarelli-Parisi ones. The notation suggests that already at the leading order (which 
is not the case for the Altarelli-Parisi kernels) their forms for timelike branchings are 
different from those for spacelike branchings. Prom ref. [12] we obtain 

X * fkSr ^ (^"- • {*}) (^"■' • ■ 

V^i I colour 

I, spin ) 

The notation 

that appears in eq. (D.l) understands that parton j © i has helicity equal to ±. 
^(ri,o) jg therefore an interference between the plus and minus helicity states of 
parton j ©i; the remaining hclicities (i.e. excluding that of parton j ©?) are summed 
over. The first term inside the curly brackets on the r.h.s. of eq. (D.l) is the ratio of 
two spinor products. We follow here the notation and conventions of Mangano and 
Parke, ref. [39]. We note that this term is a pure phase, and is therefore numerically 
well defined also in the soft limit {Ei — > 0); we have worked out its form analj^ically 
for the various phase-space parametrizations we have considered. 
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For final-state (i.e. timelike) branchings, the Q kernels read as follow [12]: 



Qgg^iz) = -ACaz{1-z), 

Qgq*(z) = 0, 
Qqq*{z) = 0. 



(D.3) 
(D.4) 
(D.5) 
(D.6) 



The corresponding results for initial-state (i.e. spacelike) branchings are: 



(D.7) 
(D.8) 

(D.9) 
(D.IO) 



Qq*g{z) = , 



Qg.q{z) = -ACf- 



The * symbol allows one to remember easily which parton is off-shell. We note that 
for branchings which involve an off-shell quark, the azimuthal term in the collinear 
limit vanishes identically. This is due to the conservation of the helicity along a 
light-quark line; see ref. [12] for details. 

We point out that in eq. (D.l) one must use the same conventions for the spinors 
in the computation of the amplitudes as in that of the pure-phase factor. As men- 
tioned above, eq. (D.l) has been worked out using the conventions of ref. [39]. On the 
other hand, in the actual numerical implementation in MadGraph the amplitudes 
are computed following the conventions of HELAS [40]. Since eqs. (D.3)-(D.6) imply 
that we need to consider only the cases in which parton j ® Hs a gluon, we can write 



Using the explicit representations of ref. [39] and of HELAS, it is a matter of algebra 
to arrive at 



where k is the gluon (light-like) four-momentum, which in the coordinate system 
chosen by HELAS has azimuthal angle equal to </7fc, and a, P are two constants, which 
are irrelevant in what follows, since the corresponding contributions to eq. (D.ll) 
vanish because of gauge invariance. Using these equations we obtain 



'°^(MP)^L"'°HmP)* = - exp (-2i (^,-ei) >lf'°^(HELAS)>l^°^ (HELAS)*. (D.14) 



We conclude this section by reporting the 4 — 2e dimensional forms of the unpo- 
larized Altarelli-Parisi kernels for z < 1 [41]: 



(D.ll) 



£''(+, k- MP) = exp {-i v9fc) £^(+, k] HELAS) + ak , 
£''(-, k] MP) = - exp (i Lpk) e^{-, k; HELAS) + /3k , 



(D.12) 
(D.13) 




(D.15) 
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Pgg{z,e) ^ Tf {z^ + {1 ^ - 2ez{l - z)) , (D.16) 
P,,{z, e) = Cf (^li^ -e(l-z)y (D.17) 

P,,{z,e) = (^i±llzii! - . (D.18) 



According to eq. (4.44), the coefficients of the 0{e^) and 0{e^) terms in these equa- 

1°' p'j,> 



tions are the kernels Pj? and Pj\^ respectively, which are used in several equations 



of this paper. 



E. Possible variants in the implementation 

As we have discussed in the text, there are a large number of arbitrary parameters 
which enter our implementation of the FKS formalism into MadGraph, and any 
physical observable is independent of them. Each parameter choice can be considered 
as a different way to implementing the subtraction procedure in a computer code. In 
this appendix we report two variants of the implementation which are not parametric 
in nature, that we may want to consider in the future, and which aim at further 
reducing the number of independent contributions to the physical cross section. We 
also present a technique alternative to that discussed in sect. 6.2 to integrate together 
the n- and {n + l)-body contributions, which is best suited for NLO computations 
matched with parton showers. 

We start by considering the contributions to the (n + l)-body cross sections due 
to {i, 1) and (i, 2) (for a given i, and assuming that both of these two pairs belong to 
^FKs)- Using the properties of the S functions, it is a matter of simple algebra to 
show that 

(B) 

X ((1 - vMM^-^'^'^) {Sa + 5,2) ^ di4yidVi<Wi . (E.l) 

Here we have set 

Vi = Vii =^ Vi^ -yi2 , (E.2) 

which follows from the fact that we work in the cm. frame of the incoming partons, 
and thus we can use eq. (4.25). Eq. (E.l) is the form used in refs. [12, 13], and has 
the virtue of allowing one to subtract the two initial-state collinear singularities in 
one single integration. On the other hand, eq. (4.29) is closer to what MadGraph 
does when it integrates tree-level matrix elements. We postpone the comparison 
between these two possibilities to a forthcoming publication. We point out that in 
refs. [12, 13] the notation was used Sf*^ — Sn + Si2, and S^j^ = Sij denoted the S 
functions relevant to both i and j in the final state. 
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We now turn to considering the contributions to the (n+ l)-body cross sections 
due to all the pairs formed by a given FKS parton i and its massive sisters. Since 
rrij 7^ for such sisters, the matrix elements are not singular in the coUinear limits. 
This implies that the distribution in yij defined in eq. (4.34) is just a regular function 
in eq. (4.29), and therefore there is a cancellation of the factors (1 — yij) between the 
numerator and the denominator. We thus immediately obtain 



j=nl '+3 
(ijOGPpKS 

E (E.4) 

(j,j)e'PFKS 



On the r.h.s. of eq. (E.3) we have changed variables yij — > yi for all j. The variable 
yi can be equal to yij for a given j, or can be defined as in eq. (E.2). The specific 
choice of yi should not matter much in terms of numerical performance, given that 
the coUinear regions yij — 1 are not associated with any peaks in the matrix elements 
(this may not be the case for particles with small masses, such as b quarks). Note 
also that multi-channel integration can be performed in eq. (E.3) following exactly 
the procedure that leads to eq. (6.29). 

In the case of processes with many massive strongly- interacting particles, eq. (E.3) 
may be a more convenient alternative to eq. (4.29), since it will reduce the number 
of independent integrations. We did not explore this possibility in the present paper, 
and we plan to do so in the future. 

We finally reconsider the problem of the simultaneous integration of the (n+ 1)- 
and n-body contributions. In sect. 6.2 we have shown how to split the n-body 
contribution into a sum of terms, each of which is then associated with a real-emission 
contribution due to a given partonic process and a given FKS pairs. We did so by 
unambiguously (but arbitrarily) "constructing" the Born-level processes by removing 
from the real processes the only gluon (in the cases in which there is a gluon) that 
is an FKS parton in Vfks- This procedure is clearly inspired by the connection 
between an (n + l)-body process and its soft limits. 

One can turn this logic around, and construct real-emission processes starting 
from a given Born process. In doing so, it is apparent that this construction can 
be achieved by considering a — > 6c branchings, with a a particle entering the Born 
process. These branchings are reminiscent of an underlying coUinear kinematics, but 
this holds only at a formal level. For example, if one has a top quark in the Born 
process, by considering the branching t ^ tg ones gets a new final state with an 
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extra gluon, which can be interpreted as a real-emission process, in spite of the fact 
that there is no colhnear singularity associated with the tg pair. 

At this point, one observes that for a given real-emission process, pairs of particles 
that lead to an underlying Born structure are in fact FKS pairs. For a given G TZn 
and a given tr G "^n+i, after constructing 7'FKs(^ii) we introduce a generalized 
Kronecker symbol 

1 it tb = r I 



V^'"^ )-U if r^^r^'V ^ ^ ^ 



We remind the reader that two processes are considered identical if one of them 
can be obtained from the other with a permutation of the identities of final-state 
particles. It is then a matter of algebra to show that 

(E.6) 

This equation suggests a definition analogue to that given in eq. (6.13). There, the 
n-body cross section rfo"*-"-* was given an (n + l)-body process as an argument. In 
this case, we give the (n -|- l)-body cross section da^"^^^^ an n-body process as an 
argument: 

da(-^^\rB)^ E E ^{rByr')4r\rn)dair\n). (E.7) 

It is quite obvious that the same procedure can be applied to the degenerate (n + l)- 
body contributions of eqs. (4.41) and (4.42), also in view of the fact that they have 
a structure identical to that of the initial-state colhnear limits of the {n + l)-body 
cross sections. We thus define: 

dar'\rB)= E E 5 {rsyf'') dair\rn) , (E.8) 
dar'\rB)= E E 5{rs,rr^)da^r\rn). (E.9) 

rR^T^n+i (i,2)G7'FKs(rH) 

da^-+'\rB) = c^ai^+'Vs) + c^^t^^+^Vb) • (E-IO) 
Equations (E.7)-(E.10) allow us to use eq. (6.15) with the formal replacement: 

E - E ■ (E-ii) 

reTZn+i relZn 
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Equation (E.7) (or, which is equivalent for the sake of this discussions, eqs. (E.8) 
and (E.9)) may appear more comphcated than its analogue, eq. (6.13). However, the 
reader should note that, owing to the presence of the Kronecker symbol, the sum over 
the {n + l)-body processes in eq. (E.7) is almost trivial, and it amounts to simply 
counting the a ^ be branchings allowed by particle identities. Given that one will 
have in any case to construct the sets of FKS pairs in order to be able to perform 
the computations of real-emission contributions, eq. (E.7) is simple to implement 
algorithmically - in particular, we stress that the current version of MadFKS does 
have all necessary ingredients to implement it. 

We conclude this session by pointing out a couple of issues. Firstly, the use of 
eqs. (E.7)-(E.10) in place of those of sect. 6.2 implies that the n-body contributions 
are computed a smaller number of times. This may lead to a significant difference 
in CPU time when including the finite term of one-loop origin V^"'^^ into MadFKS, 
whose computation for large-multiplicity final states is quite laborious. Obviously, 
the procedure of sect. 6.2 can be improved in this respect (for example, by not calling 
the function that returns V^j]J every time the Born is computed), but the technique 
described in this section is perhaps a simpler alternative. Secondly, eqs. (E.7)- 
(E.IO) are required for the matching with parton showers according to the POWHEG 
formalism. Although not strictly necessary, the same technique can also be used in 
the context of the generation of S events (which are Born-like) in MC@NLO [14]. 

F. Integration of matrix elements of given helicities 

For large-multiplicity final states, the sum over all helicities is a time-consuming 
operation. A successful strategy at the tree level is that of replacing such an exact 
sum with one performed with Monte Carlo methods, which in turn requires the 
matrix elements be available for any possible fixed-helicity configuration. This is 
indeed the case in MadCraph, and the aim of this appendix is that of showing that 
the FKS formahsm presented in this paper needs only trivial modifications to handle 
this situation. 

We recall that we compute the physical cross section using eq. (6.15), and let 
us again consider for the moment the case of e+e" collisions only, which implies 
^^(n+i) ^ g_ rpj^p _^ l)-body contribution is given in eqs. (4.28) and (4.29), while 
the n-body contribution is the sum (sec eq. (4.3)) of the quantities given in eqs. (4.4), 
(4.5), (4.12), and (4.14). In these equations, the matrix elements appear that are 
defined in eqs. (3.22)-(3.25) (where in the latter one only the finite contribution is 
actually used, as specified in eq. (B.2)). 

We now state the following fact: the equations we have referred to in the previous 
paragraph apply without any modifications to the case of a partonic process with fixed 
hehcities, except for the fact that in eqs. (3.22)-(3.25) one has to make the formal 
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replacement 

colour colour 
spin 

In other words, the sum over hehcities is now not understood in the matrix element. 

The (n + l)-body contribution has the subtracted structure given in eq. (4.38). 
In eqs. (5.17) and (5.18) we have given the necessary ingredients to construct the soft 
and coUinear counterterms respectively; these formulae however apply to the case in 
which the matrix elements are summed over helicities. When one considers the soft 
limit of the {n + l)-body matrix element for a given helicity configuration, eq. (5.17) 
is basically unchanged. The colour-linked Born matrix elements which appear on the 
r.h.s. of that equation are computed with the same helicities as those entering the 
(n-|- l)-body matrix element on the l.h.s., except for that of parton i, which is simply 
not present in the reduced process, as suggested by the notation rV The eikonal 
pre-factors are also unchanged, but need be multiplied by a factor 1/2 in the case 
in which the helicity of parton i is kept fixed; if, on the other hand, such helicity is 
summed over, no modifications are needed in eq. (5.17). 

The case of the coUinear limit is shghtly more complicated. Indeed, in the 
reduced process r-'®*'^ the helicity of the branching parton can in general assume 
any value, at fixed helicities hi and hj of partons i and j respectively. Equation (5.18) 
gets modified as follows 



+ l^^^T^^C-vMi<)M'"" (r; . (F.2) 



In the first term on the r.h.s. of this equation, h is the helicity of parton j (Bi and, 
consistently with eq. (D.2), we have denoted 

^ ' ^ ' colour 

The quantity has the same definition as in eq. (D.l), except for the sum 

over helicities, which is not performed here in accordance to the general rule of this 
section, given in eq. (F.l). We point out that in the second term on the r.h.s. of 
eq. (F.2) the sum over the helicity h of parton j (Bi is not performed, since that term 
results from the interference between the h = + and h = — states. This is also the 
reason why the kernels W do not depend on h. We finally stress that all helicities 
other than those of partons i and j have the same values in the (n + 1)- and n-body 
processes that enter the matrix elements on the two sides of eq. (F.2) respectively. 

The kernels T and W introduced in eq. (F.2) are the fully polarized versions 
of the (four-dimensional) AltareUi-Parisi and Q kernels respectively. They can be 
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computed using the results reported in app. B of ref. [12]. In the branching 
we use the following momentum fraction and helicity assignments 

a{l;ha) — > b{z; hb) c{l - z; he) , 

and wc obtain what follows (in four dimensions). 
•9 — '99 



where he — —he- 
•9 — 'QQ 



Q9 



9Q 



1 













= 0, 




= 0, 




= 0. 



l-z 

,2 



Z 



z 



q-*9q ^ ) ~ q-^qg ^ ' ' 

wii%,{z) = . 
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The T kernels relevant to the case of the branching parton with negative helicity are 

r-J:£'iz)^T:^t{z) (F.22) 

for all types of branchings; we have again used the notation hb — —hb and he — —he- 
The AltareUi-Parisi kernels at fixed helicities in four dimensions can be easily ob- 
tained from the T kernels: 

he 

Furthermore, we have 

E^a-Jc'l^) = P^^(') ' ior K = ±, (F.24) 

hbhc 
hbhc 

It is then immediate to see that eqs. (F.5)-(F.20) lead to eqs. (D.15)-(D.18) (with 
e = there), to eqs. (D.3)-(D.6), and that one recovers eq. (5.18) starting from 
eq. (F.2). In order to see this, one has to make use of the identities that relate 
AltareUi-Parisi kernels of like and unlike helicities 

= Pafi^. = PlK^. e) , (F.26) 
P;+{z, e) = P+-{z, e) ^ Ptl{z, e) , (F.27) 

which, as the notation suggests, hold in general in 4 — 2e dimensions. 

Although not relevant to the applications presented in this paper, we conclude 
this section by discussing the generalizations of eqs. (4.41) and (4.42) (the degenerate 
(n -I- l)-body contributions) to the case of processes with given helicities. These 
generahzations can be obtained from ref. [34], where the FKS subtraction method 
was extended to the case of polarized coUisions. We have 



(S) 

_^(n,o) i^^ie^. ^1 _ J_^d4>^(^\ - hy^, , (F.28) 



Afir) 



(B) 



X {r'^^'\ (1 - ^0^2) j^d<Pn{h, (1 - ^i)k2)d^i . (F.29) 
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The reduced matrix elements that appear in eqs. (F.28) and (F.29) are defined ac- 
cording to eq. (F.3); as the notation suggests, h is the hehcity of parton 1 (B i (or 
2 © i) . We are not aware of any results for the kernels K for polarized scattering, 
which are non zero in PDF schemes different from MS; therefore, when considering 
fixed hehcities we shall limit ourselves to using MS PDFs. 

The damped Altarelli-Parisi kernels at fixed hehcities, P^ab^"^^, with k — 0,1, 
have the same definition as in eq. (4.44): 

pT''' - {1 - z)P^^'^'\z) , (F.30) 

with 

pT\z) + eP^^'^iz) + 0{e') = P„V(z, e) , (F.31) 

pT%) + ^^if '(^) + 0{e') = PiJ(^, 6) , (F.32) 

where the quantities on the r.h.s. of these equations have been introduced in eqs. (F.26) 
and (F.27), and we again denoted h = —h. The 0(e°) term pj^ can be obtained 
from eq. (F.23). That equation, however, cannot be used to obtain P^^f,^"^'', since 
the kernels T have been worked out only in four dimensions. To proceed, we make 
use of the standard definitions: 

Plli^^ ^)-\ (Pabiz, e) + APabiz, e)) , (F.33) 
Piliz, ^) = l (Pabiz, e) - AP,,{z, e)) . (F.34) 

The quantities APab are conveniently introduced in the study of collisions with in- 
coming polarized beams. From refs. [41, 42] we obtain 



(F.35) 



APgg{z, e) = 2Ca (^y^ - 2^ + 1 + 2e (1 - z)^ , 

APgg{z, e) = (2^ - 1 - 2e (1 - z)) , (F.36) 

AP,,{z, e) = Cf (^i^ -e{l-z)y (F.37) 

APgg{z, e) ^CF{2-z + 2e{l- z)) , (F.38) 

after taking into account a finite scheme transformation which affects the 0{e) term 
of the qq kernel [42] . 
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